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Introduction 

The structural dynamicist’s areas of responsibility require interaction with most other 
members of the wind turbine project team. These responsibilities are to predict structural 
loads and deflections that will occur over the lifetime of the machine, ensure favorable 
dynamic responses through appropriate design and operational procedures, evaluate potential 
design improvements for their impact on dynamic loads and stability, and correlate load and 
control test data with design predictions. Load prediction has been a major concern in wind 
turbine designs to date, and it is perhaps the single most important task faced by the 
structural dynamics engineer. However, even if we were able to predict all loads perfectly, 
this in itself would not lead to an economic system. Reduction of dynamic loads, not 
merely a “design to loads” policy, is required to achieve a cost-effective design. 
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The two processes of load prediction and structural design are highly interactive: loads 
and deflections must be known before designers and stress analysts can perform structural 
sizing, which in turn influences the loads through changes in stiffness and mass. Structural 
design identifies “hot spots” (local areas of high stress) that would benefit most from 
dynamic load alleviation. Convergence of this cycle leads to a turbine structure that is 
neither under-designed (which may result in structural failure), nor over-designed (which 
will lead to excessive weight and cost). 

This chapter introduces some of the physical principles and basic analytical tools 
needed for the structural dynamic analysis of a horizontal-axis or vertical-axis wind turbine 
(HAWT or VAWT). This is done through discussions of two subjects that are fundamental 
building blocks of our understanding of the structural-dynamic behavior of wind turbines: 

— Single-degree-of-freedom dynamic load model of a HAWT blade, following the 
development of the FLAP computer code at the National Renewable Energy 
Laboratory [Thresher et al. 1986, Wright et al. 1988]; 

— Theoretical and experimental analysis of the vibration modes of a VAWT 
system, following methods developed at the Sandia National Laboratories 
[Cam tetal. 1982], 

These discussions are written for engineers familiar with structural mechanics, but no 
specific knowledge of wind turbine rotors is required. Other technical disciplines with 
which the wind turbine structural dynamicist should become familiar are aerodynamics 
(required to determine wind forces and aeroelastic effects), rotor dynamics (with its own 
special set of accelerations and responses), statistics (required in dealing with wind 
turbulence and fatigue life prediction), mathematical modeling, and testing. 

Role of Structural Dynamics in the Overall System Design 

The design cycle in a typical wind turbine project is usually divided into concept 
definition and detail design stages. During the concept definition period trade studies are 
conducted to determine the overall configuration. The structural dynamics engineer is 
concerned most with the relative load levels and the operational and technological risks 
associated with the various design alternatives. Initially, only the feasibility of each concept 
is judged, with little time spent on refinements. For example, aeroelastic stability might be 
addressed only for unconventional concepts where experience is lacking. 

In the detail design stage the structural dynamicist must furnish the design loads and 
determine design and operational requirements that insure acceptable dynamic behavior. 
These requirements might include drive train damping and spring rates, balance tolerances, 
stiffness ranges for bearings, maximum operating wind speed, allowable yaw misalignment, 
yaw rate limits, etc., all based on the results of analysis using various mathematical models. 

The key objectives that the design team seeks to achieve by relying on the guidance 
and skill of the structural dynamicist include the following: 

— To select a configuration and design approach which alleviate dynamic loads 

— To minimize the sensitivity of dynamic responses to inevitable differences between the 
“as-designed” and “as-built” physical properties of the structure 

— To place natural frequencies favorably with respect to turbine operating speeds 

— To insure aeroelastic stability, without which safe operation is impossible 

Accomplishment of the first of these objectives must be based on the team’s collective 
understanding of the dynamic behavior of various wind turbine configurations. It is guided 
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by past experience and frequently involves trade-offs with factors outside the realm of 
structural dynamics ( e.g ., cost, manufacturing, and energy capture). The remaining key 
objectives are accomplished by straightforward (though sometimes complex) mathematical 
modelling of the selected configuration. Of course, the likelihood of a good detail design 
solution should be evident at the completion of the concept definition stage. 

Wind Turbine Substructures and Subsystems 

General background material on wind turbine nomenclature, configurations, and major 
substructures may be found in Chapters 2 and 3. It is its rotor which sets a wind turbine 
apart from other structural systems. Therefore, the thrust of this chapter is the understand- 
ing and analysis of wind turbines rotors. The structural dynamic behavior of the power 
train and support structures of a wind turbine can usually be analyzed by well-established 
methods, so these subsystems will not be discussed here. 

Mathematical Models 

While certain analyses can be accomplished with simple formulas, experience has shown 
that the analyst will need separate, specialized computer models to determine 

— System structural modes and natural frequencies which incorporate rotational effects 

— System loads, both steady-state and transient, incorporating all important degrees of 
freedom; includes deflections, vibrations, and power output quality 

— Aeroelastic stability, including the effects of structure/control system interaction. 

In theory, one comprehensive computer code could be used for all three of these tasks, but 
it might not be as efficient to develop or to run as separate specialized codes. 

Modal Analysis Models 

Modal analysis is the determination of the set of discrete patterns or mode shapes, each 
with at its own modal frequency and modal damping, which a vibrating structure describes. 
These patterns are also known as natural modes. Modal analysis of a wind turbine is 
generally based on well-developed finite element procedures (e.g. NASTRAN). Each 
substructure in the turbine system is modeled separately to determine its own mode shapes 
and frequencies. With special care, substructures can be coupled together to produce system 
modes. Modal models can be derived directly or spawned from more-detailed stress 
analysis models using standard modal reduction techniques. For rotors, centrifugal 
stiffening and gyroscopic effects can be incorporated directly into the finite element code, 
or they can be computed externally and applied later. 

System Loads Model 

The system loads code should be tailored to suit a specific wind turbine configuration, 
since a code that attempts to analyze all possible wind turbine types can become unwieldy. 
For example, the important dynamic loads usually result from motion in only a few 
dominant modes that are characteristic of a particular configuration, so the degrees of 
freedom used in the system loads model are generally a selected subset of the natural modes 
computed with the finite element models described above. While it is possible to use all 
of the finite-element degrees of freedom in the loads model, a modal approach is commonly 
used for analyzing system loads. 
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Direct numerical integration of the forced equations of motion in the time domain (i.e., 
output is load vs. time) is the most straightforward and informative method of solution. 
This procedure, which is often referred to as system simulation, handles both steady-state 
and transient response. Since the governing equations of motion contain time-varying 
coefficients and are generally nonlinear, frequency domain techniques have limited value. 
Furthermore, seeking closed-form analytical solutions to these complex equations is usually 
unnecessary with the capabilities of today’s digital computers. The engineer’s time is best 
spent in understanding the results of a system simulation and then exploring methods for 
reducing dynamic loads. 

Aeroelastic Stability Models 

Aeroelastic stability analysis is often separate and distinct from the modal and loads 
analyses for two reasons: First, it is preferable to linearize the equations of motion and 
examine the resulting eigenvalues for stability, rather than to examine a myriad of time 
histories to determine if the response is stable. Of course, if the mechanism of instability 
is highly nonlinear, simulation in the time domain may be the only recourse. Most potential 
wind turbine instabilities can be adequately assessed using linearized models. 

The second reason for a separate model is that the modes which are important for 
aeroelastic stability often differ from those needed for loads analysis. In a successful 
system, instabilities will occur outside the limits of planned operation. Hence, a blade 
torsional mode which incites instability at a high rotor speed may well be dormant during nor- 
mal operation and, hence, be omitted in the loads analysis. It should be noted that aeroelastic 
stability analysis, like modal analysis, can be carried out with a finite-element model. 


Dynamic Load Model of a HAWT Blade 

A beam model of a wind turbine blade is generally suitable for structural-dynamic 
analysis. It will differ from the small-deflection theory beam models used in conventional 
analyses of non-rotating structures, however, because the axial loads on the blade 
significantly effect the lateral and torsional deflections. In this respect it is more like the 
beam-columns representations used in elastic stability analysis. General three-dimensional 
theory is quite complicated and has provided fodder for more than one doctoral thesis. It 
is still subject to controversy as to which terms are important and which are not. Rather 
than beginning with the general case, we will develop the equations for uncoupled flapwise 
{i.e. out-of-plane or flatwise) motions to acquaint the reader with the physics involved. In 
many cases, these will be suitable for the analysis of HAWT rotor loads. 

Elastic blade flapping equations appear in many sources and are in a sense, classical. 
A relatively brief derivation is given here to highlight the approximations and assumptions 
that are implicit in the equations which are commonly used. We begin by assuming that 
the strains everywhere in the blade are small (less than 10%) and below the elastic limit. 
Other assumptions will be introduced at appropriate points in the derivation. 1 Equations of 
motion for wind turbine components are generally derived using a “strength of materials” 
rather than a “theory of elasticity” approach. 


1 Note here that the assumptions of small strain and elasticity do not mean that the 
deflections must be small. Everyone is aware that a steel ruler can be bent into an arc with 
deflections much greater than its cross-sectional dimensions, and yet it will spring back 
elastically to its original shape. 
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Coordinate System Definitions 

Figure 11-1 shows the orientation of the HAWT blade under analysis with all the 
intermediate coordinates required to represent the blade motion. The X-Y-Z coordinates are 
the fixed reference system. The mean wind velocity at the hub elevation, V hub , and its 
fluctuating components — SV X , &V y , 5V Z — are given in this system. The rotor axis may 
be tilted at a fixed angle %■ It may also have a prescribed, time-dependent yawing motion 
given by the yaw angle <()(/). The yaw axis is coincident with the Z-coordinate axis. The 
center of the rotor hub, H, is located a distance a from the yaw axis. The blade may be 
rigid from the axis outward to a distance h, which may also be interpreted as the outer 
radius of the hub. The blade airfoil shape may begin at h or at a station further out along 
the blade z axis. The blade may be coned at a prescribed angle (5 0 , as shown in the figure. 

The x,y,z coordinates are located in the surface of revolution that a rigid blade would 
trace in space, with the y-axis normal to this surface. The x p -y p -z system are the blade’s 
principal bending coordinates, where the z p -axis is coincident with the elastic axis of the 
undeformed blade. Bending takes place about the x p -axis. It is further assumed that the 
blade principal axes of area inertia remain parallel along the z p -axis. Any influence of 
blade twist on bending displacements is neglected. As a reference for computations, the 
angle 0 p is set equal to the orientation of the principal bending plane relative to the cone 
of rotation (i.e., the x-y-z system) at a selected station of primary interest along the blade 
span. The third coordinate system illustrated in Figure 11-1 is the system, which 
is on the principal axes of the deformed blade along each point along the elastic axis. 



Figure 11-1. Fixed, rotating, and deformed-blade coordinate systems. Positive 
displacements and rotations are shown. [Wright et al. 1988] 
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The coordinate transformation that takes vector components from the fixed X-Y-Z 
system to the rotating x p -y p -z p coordinates of the undeformed blade is as follows: 


x p 

y P 

, Zp . 
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( 11 - 1 ) 


This transformation has been linearized by assuming a small yaw angle <|>, a small tilt 
angle %, and a small coning angle P 0 . The angle of the principal axis, 0 p , may be large. 
Of course, the blade azimuth angle, \|f, can vary from 0 to 27t radians. The inverse trans- 
formation taking the x p -y p -z p components into the X-Y-Z system is given by the transpose 
of the 3x3 matrix in Equation (1 1-1). The blade position for \|f = 0 is straight up. 

Assuming small deformations, the transformation for obtaining rj-^-2, components is 
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( 11 - 2 ) 
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The necessary intermediate coordinate transformations for components in the x-y-z, 
x-y-z r , and x h -y h -z h coordinate systems shown in Figure 11-1 can be obtained as follows: 

— Set 0 p = 0 for transforming to the x-y-z system; 

— Set 0 p = |3 0 = 0 for transforming to the x r -y r -z r system; 

— Set 0 p = P 0 = V|f = 0 for transforming to the x h -y h -z h system. 


Moment-Curvature Relationship 

The HAWT blade is assumed to be a long, slender beam so that the normal strength-of- 
materials assumptions concerning the bending deformations are valid. Figure 1 1 -2 shows 
an infinitesimal element of the deformed blade. It is assumed that the blade bends only 
about its weakest principal axis of inertia, which is about the x p -axis in the figure. No other 
deformations are considered in this analysis. The strength-of-materials model for elastic 
bending assumes a one-dimensional form for Hooke’s Law that neglects all stresses except 
the longitudinal (in this case spanwise) bending stress. It results in the following 
differential equation relating bending moment to curvature in our blade: 


M = - El < ^ V - 
x ' p x ' p dzl 


- EI x<p v 


(11-3) 


where 

M = bending moment about the blade principal axis of inertia, x (N-m) 
[ ]' = d[]/dz p , etc. 

E = modulus of elasticity (N/m 2 ) 

I = area moment of inertia of the airfoil section about the x -axis (m 4 ) 

x,p p V 7 

v = displacement in the y (flap) direction (m) 



Figure 11-2. Deformed blade element showing forces and moments, all acting in a 
positive sense. [Wright et al. 1988] 
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Equilibrium Equations 

The equations of equilibrium are derived by summing forces and moments in the three 
coordinate-axis directions and equating the sums to zero. Referring to Figure 11-2, 
summing of forces and moments gives the following six equations: 


Vx,p + Px,p = 0 

(ll-4a) 

Vy,p + Py<p = 0 

(11 -4b) 

^ + Pz,p = 0 

(11 -4c) 

p - Vy.p + Tv ' + Qx,p = o 

(ll-5a) 

My,p + Vx,p + Qy,p = 0 

(11 -5b) 

Kp - V *,p v ' + 9z, P = 0 

(11 -5c) 


where 

p x p , p y p , p zp = applied (external) force loading per unit length (N/m) 
q x , q y p ,q zp = applied (external) moment loading per unit length (N-m/m) 

M ,M ,M = internal bending moment loads (N-m) 

V ,V , T = internal force loads (N) 

*>p y>p N 7 

The terms load and loading are often used interchangeably in the literature, so some 
clarification of their use in this chapter (and generally in this book) is in order. We shall 
use the term loading to describe an external action applied onto the structure. Thus, p and 
q are the sum of all aerodynamic pressures and inertial body forces and are here called 
loadings. The internal moment and force responses of the structure to these loadings - M, 
V, and T - are termed loads. It has become the general practice in the structural analysis 
of wind turbines that if the term “load” stands alone it refers to internal moment and force 
responses. 

Differentiation twice of Equation (11-3) and once each of Equations (ll-5a) and (11- 
5b) allows these three equations to be combined with each other and with Equations (11 -4a) 
and (11 -4b) to eliminate M xp , V xp , and V . In addition, Equations (1 1 -5b) and (ll-5c) can 
be combined to eliminate the shear force load V . These eliminations reduce the previous 
seven equations to the following four combined moment-curvature-equilibrium equations 

in the four unknowns v,M ,M , and T: 

x,p y.p 


Flapwise Bending : 

( — v "EI yp )" + (T v'Y + q' xp + Py , p = 0 

(ll-6a) 

Edgewise Bending : 

^ y.p + Qy,p ~ Px,p = 0 

(11 -6b) 

Pitchwise Torsion : 

M ZtP + MypV + qy,pV + q Zt p = 0 

(ll-6c) 

Spanwise Tension: 

T + Pz,p = 0 

(11 -6d) 
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Wind Speed Models 


Wind Shear Gradient 


The vertical gradient in steady wind speed with elevation above and below the elevation 
of the rotor axis is often modeled by a simple power law (see Chapter 8) as 


U(z) = U„(zlHy 


(11-7) 


where 


U(z) = free-stream wind speed at elevation z (m/s) 

U H = free-stream wind speed at hub elevation (m/s) 
z = elevation above ground level (m) 

H = elevation of the rotor hub (m) 

m = empirical wind shear exponent (same as a in Eq. (8-1 1)) 


The wind speed gradient may be described in polar coordinates centered at the hub 
elevation by a binomial series, as follows: 


U(z) = U„ ( rCOS ^ + g ) = U„[ 1 + W s (r,y)] 

, . ( r\ m(m - 1) / r \ 2 2 

Ws(r, \|f) ~ m cost)/ + cos z \|t 


( 11 - 8 ) 


(11-9) 


where 


r = radial distance from the rotor axis (m) 
W = wind shear shape function 


In this series expression, powers of r/H higher than 2 have been neglected. 

Tower Shadow Deficit 

The presence of the tower causes a deficit in the wind speed, with the downwind 
reduction often referred to as tower shadow and the much-smaller upwind reduction as a 
bow wave effect. Here we shall refer to both as tower shadow. Neglecting any effects of 
retardation by the rotor, we will assume that the deficit has a spatial distribution of the form 

5Vr(z,V) = -W r (\|0t/(z) 


to + tp COS[p(\|f - 7t)] 

WKy) = for n — Vo < V|r < Jt + \|fo 

0 elsewhere 


( 11 - 10 ) 


p = &7t/\|ro 


( 11 - 11 ) 
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where 

8V t = wind speed deficit caused by the tower (m/s) 

W T = tower shadow shape function 
tg, t = empirical scale constants 
t|t 0 = half-angle of tower shadow sector (rad) 
k = number of waves in the tangential profile of the shadow 

The parameters t 0 , t p , \|t 0 , and k are selected to give the desired approximation for the 
tangential profile of wind speed in a pie-shaped region with a central angle of 2\|/ 0 , centered 
on a blade azimuth angle of ji. For example, the shadow of a shell tower is often modeled 
as a sin 2 function using the following parameter values: 

t 0 = t = one-half the maximum ratio of deficit to free-stream wind speed 
k = 1 

A truss tower with three legs produces a shadow with three peaks, and this is modeled by 
selecting k = 3. 

Spatial Turbulence Model 

As discussed in Chapter 8 and illustrated in Figures 8-16 and 8-29, small-scale 
turbulence (i.e. having significant spatial variations in wind speed within the swept area of 
the rotor) will cause a moving rotor blade to experience a wind power spectrum with peaks 
at multiples of the rotor speed. In effect, the blade will be subject to harmonic forcing 
functions with frequencies at integer multiples of its rotational speed, and it will respond 
dynamically to them. Estimating the amplitudes of these turbulence forcing functions is 
beyond the scope of this chapter and is, in fact, the subject of current research on wind 
characteristics. Here we will make provision for wind turbulence excitations and emphasize 
that the harmonic content of these excitations plays an important role in determining the 
structural-dynamic response of wind turbine blades. 

We shall assume that the distribution of turbulence-induced variations in wind speed 
across the rotor swept area changes slowly compared to the period of rotor rotation, so our 
model will be non-uniform in space but uniform in time. In the general case illustrated in 
Figure 11-1, the wind turbulence will have components in all three of the X-Y-Z directions, 
and each of these may vary with position in the rotor disk area. Thus 


5t/(r, \|f, t) = 


8u(r, \|t, 0) 
8 v(r, \|t, 0) 
8w(r, \|t, 0) 


( 11 - 12 ) 


where 

8U(r,\\t,t) = turbulence in the free-stream horizontal wind (m/s) 
5i<(r,\|/,0) = spatial turbulence in the X direction (m/s) 

8v(r,\|/,0) = spatial turbulence in the Y direction (m/s) 

8w(r,\|/,0) = spatial turbulence in the Z direction (m/s) 
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Induced (Retardation) Wind Speeds 

In an unducted rotor, extraction of power reduces the speed of the wind at the rotor 
blades compared to the free-stream wind speed, as discussed in Chapter 5. This reduction 
is termed the induced wind speed and its magnitude is usually given as a fraction of the 
free-stream wind speed by a parameter called the axial induction factor (see Fig. 5-11). 
Calculation of axial induction factors is routinely performed using aerodynamic performance 
models, balancing the rate of change of the momentum in the wind stream with the rotor 
thrust. Because thrust is not normally uniform across the swept area of the rotor, advanced 
performance codes determine radial and circumferential distributions of the induction 
factors. 

Some dynamic-load computer models also include the calculation of induced wind 
speeds. Whatever their source ( i.e ., external or internal to the load model), provision must 
be made for induced wind speeds, because they can have a significant effect on airfoil 
angles of attack and on aerodynamic loading. In our wind model this is done as follows: 

Vi(r,y)= a(r,y)U H (H-13) 

where a(r,i) = prescribed spatial distribution of axial induction factors 
Combined Wind Effects 


The net wind speed at point in the plane of rotation is the sum of free-stream, wind 
shear, tower shadow, turbulence, and induction components. These are usually combined 
in a linear fashion, neglecting minor interaction terms. To help avoid ill-conditioned 
simultaneous equations later, wind speeds are then normalized by the tip speed of the rotor. 
The net wind velocity vector is specified in the X-Y-Z coordinate system and is composed 
of steady and turbulence terms, as follows: 



0 


8h 

V W (X- Y - Z) = RC2 

Vr + 6V, 

+ iRQ 

5v 


0 


_8h». 


(ll-14a) 


V r 


V H [l - a(r, t|Q] 
RQ. 


(ll-14b) 


8V, 


VH\Ws(r, y) - W T (r, V|Q] 
RCl 


5 u = 8 u(r,\\r)/RCl 
5v = 5v(r, \|f)//?S2 
8»v = 5w(r,\|f)//?Q 


(ll-14c) 

(ll-14d) 
(1 1 -1 4e) 
(ll-14f) 


where 

Vw(r, t) = net wind speed at the plane of rotation (m/s) 
R = rotor tip radius (m) 

Q = rotor rotational speed (rad/s) 


ASM_Wind T urbine_Ch1 1 .indd 


617 


MTC 


02/16/2009 10:06AM 



618 


Chapter 11 

Kinematics of the Blade Motion 


Velocity Analysis 



We will now derive equations for the velocity of a point on the deformed blade in the 
r|-£-^ deformed blade system, starting in the X-Y-Z ground-based coordinate system and 
then transforming the results first into x p -y p -z p principal blade coordinates and finally into 
the coordinates. Referring to Figure 11-1, we can illustrate the required velocity 

vector analysis as follows: Designate the X-Y-Z coordinates as the a reference frame. Call 
the x p -y -z principal blade coordinates, located at a point B on the elastic axis of the blade, 
as the $ reference frame. The velocity of the same point on the deformed blade in the a 
system, here designated as A, may then be written symbolically as 

Va(oc) = 

d 

■■ V B (a) + —r m (P)+ Q xr m 

Ql 

where subscript AJB indicates a vector from point B to point A. Performing the computations 
indicated by this equation and transforming the result the result using Equation (11-1) and 
then (11-2) gives 

Va(x\ - t, - 5) = *£2 

rcQp - <j)[(d + Po r)cQ p + v]cvjr - <t>(r£0p),n|/ 
rsQp + v - <j>[(d + Po»*)50p]c\|r + <j) (rc0 p )s\jr 
(v'r - v)sQ p - <j>[rf + (v - v'r)c0 p ] s:\jr 

(11-1 6a) 


r 

r ~ R 

(11 -16b) 



(ll-16c) 


‘-i 

(11-1 6d) 


V 

v - R 

(ll-16e) 

where <j) = yaw angle of the rotor axis (rad) 

d = distance from the tower axis to the rotor hub (m) 


and bold print designates a dimensionless variable. In addition, it has been assumed that the 
order of magnitude for the various terms are as follows: 

Order 1 variables : r 

1/2 

Order e 0 variables : <j) 

Order Eq variables : d, Po, v, v, v' , H/R 


Terms of order e 0 2 and higher have been neglected in Equation (1 l-16a). 
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The relative velocity of the wind with respect to the moving rotor blade is computed 
by transforming the wind velocity components from Equation (11-14) into the deformed 
blade coordinates and then subtracting the blade velocity of Equation (1 1-1 6a). This gives 


Vrd 

RCl 



(ll-17a) 


= -rcQ p - (V, + 5 V, + 8 V ,)sQ p + c\|t[(8V* + $V r )cQ p 

fill 

+ 0{(d + Por)c0 p + v}] 


(ll-17b) 


= -r S e p -v + ( Vr + 5 Vr + 8 V y )cQ p + c\|/[(8F* + <\ )V r )sQ p 
Ril 

+ i>{(d + $ 0 r)sQ p }] (ll-17c) 

+ 5t|t{ -(SV Z - lVr)sQ p - $(rcQ p )} 


- (vV- V)50 p + Ct|/(8V Z -%Vr) + ( Po + VcQpWr 
+ ,«|r[8V*+ ifV r + ^{rf+(»- vV)c0 p }] 


(ll-17d) 


where 

V rd = relative wind speed at airfoil section in ly^-E, coordinates (m/s) 

% = relative chordwise wind speed at airfoil section (m/s) 

Vj. = relative normal wind speed at airfoil section (m/s) 

V/ = relative spanwise wind speed at airfoil section (m/s) 

Terms of order e 0 2 have been discarded in these equations. The following assumptions 
have been made regarding the order of the velocity components: 


Order 1 velocities: V T 

Order £q velocities : 8V r , 81*, 81^,, 8V^ 


Acceleration Analysis 

Referring again to Figure 11-1, our acceleration analysis will recognize the fact that 
the mass of the blade is not all concentrated on the elastic axis but is distributed across the 
airfoil section. Acceleration equations will, therefore, be derived for an arbitrary point P 
at coordinates (r|,Q within the cross-section. The point A in the velocity analysis is the 
same as P(0,0). Using the same <x-P reference frame designations as for the velocity vector 
analysis, the acceleration of P is given by the usual five-term acceleration equation as 


ASM_Wind T urbine_Ch1 1 .indd 


619 


MTC 


02/16/2009 10:06AM 



620 


Chapter 11 


ap(cc) - a B (a) + ^-xr P/B + 2£lx^r P/B (fi)+ £lx(£2xr P/B ) + ^ r P/B (11-18) 


The indicated operations of Equation (11-18) must be carried out, and the results must then 
be transformed to the P coordinate system and linearized. This tedious activity gives the 
expressions for the acceleration components in the x p -y p -z p coordinates. For the small- 
deformation theory of this analysis, there is no difference between the x -y -z and q-£-£ 
coordinate systems with respect to structural equations. For convenience, therefore, we will 
use the same q-£-ij notation for the accelerations as for the relative velocities. This gives 


a P = a p ^ 


(ll-19a) 


ap,Ti = -5'0 p (rQ 2 Po + 2<j> Qrc\|/ + $rsv|/ + vQ 2 c9 p ) 
- t^tfcQpsQp - qQ 2 c0 2 


(11-1 9b) 


a P ^ = - cQ p (rQ 2 Po + 2 <j)Qrc\|/ + cb rvqr ) - vQ 2 c0 2 
- £ Q 2 r0 2 - q Cl 2 sQ p cQp + v 


(11-1 9c) 


a P ^ = -r£2 2 -2Qvs0 p 


(1 1-1 9d) 


where 

a pj] , a vt , a vl = components of the acceleration vector at an arbitrary point P within the 
blade, in the normal, chordwise, and spanwise directions, respectively 
(m/s 2 ) 

Aerodynamic Loading 

Referring to Figure 6-1, the lift and drag forces on an airfoil section are given by 

dL=\ 9 C L cV} el dl 


dD = V -?C D cV} el d\ 


where 


dL = increment of lift force, normal to the relative wind (N) 
dD = increment of drag force, parallel to the relative wind (N) 
p = air density (kg/m 3 ) 
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C L , C D = lift and drag force coefficients, respectively; prescribed functions of a 
a = angle of attack, from the relative wind to the chord line (rad) 
c = chord length of the airfoil section (m) 
dt, = increment of spanwise length (m) 

We see from Figure 6-1 that the directions of the lift and drag forces are determined by the 
direction of the relative wind and not by the blade geometry. As a result, lift and drag 
directions will vary along the span of the blade and change with changing wind and 
operating conditions, which is not convenient for structural analysis. 

The force components of interest for structural analysis are those related to the principal 
inertia axes of the entire blade. To locate these directions, we first define a reference 
station along the spanwise axis of the blade. The reference station is any location of special 
structural or aerodynamic interest. Referring to Figure 11-3, the angle between the cone 
of rotation and the chordline at this reference station is defined here as the pitch angle, Qp. 



Figure 11-3. Velocity triangle at the blade reference station, and coordinate systems 
used for computing aerodynamic loading. The twist angle is zero. [Thresher et al. 1986] 

The chord line at the reference station becomes the rj-axis. Next, we define a reference 
plane as the surface generated by a rj-axis line moving along the spanwise axis, remaining 
parallel to the chord line at the reference station. Finally, we define the twist angle, 0 ( , as 
the angle between the local airfoil chord line and the reference plane. The force 
components of structural interest can now be expressed by 
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dA r, = ±pcV rel [C L ( a)V c - C D (a)V n ]d^ 


(11 -20a) 


dA c = ipcV^tCiCa^-n C D (a)V c ]^ 


(11 -20b) 


a = 0^ - 6f 


(11 -20c) 


Or, = arctan(V^/V n ) 


(ll-20d) 


where 

= increments of aerodynamic loading normal to and parallel to the blade 
reference plane (N/m) 

0^ = angle of the relative wind from the blade reference plane (rad) 

0 ( = local twist angle relative to the blade reference plane; positive twist 
rotates the trailing edge downwind, as does positive pitch, B p (rad) 

Most HAWT blades and their pitch control mechanism (if they have one) are relatively stiff 
in torsion about the spanwise axis, so aerodynamic pitching moments can be neglected. 
However, if a blade has pitch flexibility (e.g. one that is self-twisting for controlling power 
at high wind speeds) pitching moments must be included in the loading analysis, as follows: 


9i ~ Vrel C m, a.c. €e.a.cdA^ 


(11-21 a) 


(ll-21b) 


increment of aerodynamic pitch moment loading (N-m) 
pitch moment coefficient about the aerodynamic center of the airfoil 
relative eccentricity of the elastic center of the airfoil; positive for 
elastic center between leading edge and aerodynamic center 
distance from leading edge to aerodynamic center; usually c/4 (m) 
distance from leading edge to elastic center (m) 


Inertia Loading 

The distributed inertia forces and moments acting along a unit length of the blade can 
be computed using Newton’s laws as follows: 


^ 1 ^ 

Pi = J -a p dm = J J - a P p b dy\d £ 


(11-23) 


?/ = J —r c x apdm = J J (r c x a P )p b dr\dC, 


(11 -24a) 
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r c = 




L o J 


(11 -24b) 


where 

dm = mass of an incremental length of blade (kg) 
p h = average mass density of the blade section (kg/m 3 ) 

If the force of gravity is considered as deriving from an acceleration in the -Z direction, it 
can be treated in a similar manner. 

Distributed Loading on the Blade 

The combined loading on a HAWT blade caused by aerodynamic forces, shaft rotation, 
yaw motion, and gravity are 


Pr\ = dAyy + m[(rQ 2 Po + 2Q^rcv|t + <j>mj/)s0p + vQ^GpcGp + e 11 Q 2 c0 2 


+ fngixsQp + c Q p .s \j/ - p o ,s-0pC\j/) 

(11 -25a) 

p ^ = dA ^ + m[-(rQ 2 p 0 + 2Q<j)rc\|f + <j>m|/)cG p + vQ 2 sG 2 
+ e n £2 2 sQ p cQ p - v] + mg(~xcQ p + s-GpS'Vj/ + P o c0 p c\j/) 

(11 -25b) 

p ^ = m(rS2 2 + 2vilsQ p ) - mgctjf 

dn = v'Q 2 sGX 

(11 -25c) 

H — mrCl 2 e r, - v'Q}cQ p sQ p I^ + mge^cy 

(ll-25d) 

= %,a + n 2 sQ P cQp(I^ - /” n ) + mgeyysyQp 

(ll-25e) 

m = J Pbdr\dC, 

Section 

(ll-25f) 

gT1 = m J" J" piT|flfr| ^ 

Section 

(ll-25g) 

J J 

Section 

(ll-25h) 

T y\x\ = J J Pbf\ 2 dr\dt; 

(ll-25i) 


Section 
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where 

m = blade mass per unit length (kg/m) 

% = ^-coordinate of the section center of mass (m) 

I^ m = minimum mass moment of inertia of section, per unit length (kg-m) 

= maximum mass moment of inertia of section, per unit length (kg-m) 

Governing Equation of Motion 

For this single-degree-of-freedom system, with flap displacement v as the only dynamic 
motion, the standard equation of motion for an incremental length of the blade has the form 

[M] v + [J]v + [AT] v = [L] (11-26) 


where 

[M] = mass function (kg/m) 

[J] = damping function (N-s/m 2 ) 

[K] = stiffening function (N/m 2 ) 

[L] = loading function (N/m) 

To develop this equation governing flap motion, it is necessary to begin with Equations (11- 
6) defining the equilibrium-moment curvature relationships for a blade element. Equation 
(ll-6a) for flapwise bending will be converted into an equation of motion. 

Spanwise Tension 

The spanwise tension, T, in Equation (11 -6a) can be obtained by directly integrating 
Equation (1 l-6d) to give 


R 

T{\) = J m(rQ 2 + 2v£ls% - gcy)dt, (11-27) 

$ 

For the integration limit, we neglect any difference between the tip radius R and the length 
of the blade. 

Modal Flap Displacement Model 

In order to reduce the flap motion equation to an ordinary differential equation required 
for a computerized solution, a modal model is used, in which the flap displacement is 
assumed to be of the form 

v£,0 = 2>*(0 Y*(5) k = 1 , 2, ... (11-28) 

k 

where 

k = number of mode shapes in the displacement model 

s k = time-dependent displacement scaling factor for the Mi mode (m) 

y k = spanwise shape function for the Mi mode 

The spanwise shape functions, which must satisfy the blade kinematic and natural (force) 
boundary conditions, are usually prescribed as the mode shapes for free vibration of the 
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blade. The time-history of sfi) is determined by numerical integration of the equation of 
motion during consecutive revolutions of the rotor, until each s k and its first time derivative 
at \|f = 2n are equal to the same parameters at \|/ = 0, within a set tolerance. This is the 
so-called trim solution. 

The flap mode shapes of importance to the structural-dynamic analysis of HAWT rotors 
can be divided into the following two classes: (1) dependent blade modes, in which the 
motions of a given blade are constrained or influenced by the motions of another blade in 
the rotor, and (2) independent blade modes, in which the motions of a given blade can be 
analyzed separately from the other blades. An example of a dependent mode is the rigid- 
body teetering motion of two-bladed HAWT rotor, in which the mode shapes of opposite 
blades must always be antisymmetric and without curvature at the hub. Dependent blade 
motions will not be discussed further here because of the additional complexity introduced 
by the constraints at the hub. The reader is referred to Wright et al. [1988] and [Wright 
and Butterfield 1991] for the equations required to perform a dynamic analysis of a teetered 
HAWT rotor and sample comparisons with test data. 

The two classic examples of independent blade modes are the cantilevered shape, in 
which both v and V are zero at the hub, but v" may be non-zero; and th e flapping hinge 
shape, in which both v and v" are zero at the hub, but v' may be non-zero [Spera 1975], 
In the remainder of this analysis, we will limit ourselves to one independent flapwise mode, 
for simplicity (i.e., k = 1). Fortunately, experience has shown that one cantilever flapwise 
mode is usually sufficient for the calculation of blade loads in a rigid-hub rotor. Moreover, 
blade loads are much more sensitive to the scale of the wind input models (for wind shear, 
turbulence, and tower-shadow) than they are to higher blade modes. The exception to this 
general observation is the case in which the natural frequency of the next-higher mode is 
equal to an integer multiple of the rotor speed, and a resonance condition may be present. 

Equation of Motion for One Independent Blade Mode 

Substitution of the assumed form of the blade displacement from Equation (1 1-28) into 
Equation (1 l -6a) with k set equal to 1 gives 

[- (syfEI^f + (TsY)' + <4 + = 0 (11-29) 

Now, using a Galerkin approach [e.g. Dym and Shames 1973] gives us 

R 

J {[- (rtf’EItf’ + (Tsy'y + q\ + pc] yd\ = 0 (1 1-30) 

o 

Integrating this equation by parts and using the boundary conditions required of y(ff) lead 
to an ordinary differential equation in s(t), which can then be integrated numerically in the 
time domain. 

The numerical integration approach is best explained by rearranging Equation (1 1-29) 
so that the acceleration term remains on the left-hand side but all other terms are moved to 
the right-hand side. If we perform our time-domain integration in time steps of duration 
dt, we can now calculate the acceleration at time t in terms of the right-hand side of the 
equation evaluated with calculations made previously at time t - dt. We continue this 
forward-integration process in time until convergence to a trim solution is obtained. As 
discussed previously, trim is achieved when the flap displacement and velocity at \|/ = 2tt 
are equal to the same parameters at \|/ = 0, within a specified tolerance. 
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Carrying out the integration of Equation (11-30) by parts gives the following 
expressions for the flap accelerations of an independent blade: 

s'M = - sKb — s (Q. 2 Kq + IClsQpsKc - c\\tgK g ) 

{Bending) {Tension Stiffening) 

-cQ p (fl 2 + 2Q<j)c\|/ + §s\\[)Mr 
{Rigid -Body Motion) 

+sCl 2 sQ 2 K q 


{Inertia Moment Stiffening) 

+ Q 2 sQ p cQ p Mb + F a 

{C.G. Imbalance) {Aero Forace) 

+ s£1 2 sQ 2 M 

{Inertia Force Stiffening) 

+ g(-%c®p + + p 0 c8 p c\|/)Af g 

(Gravity) 


(11-31) 


The various coefficients in this acceleration equation are given by the integrals in Equations 
(11-32) which are evaluated numerically. A trim solution of the equation of motion is then 
found by forward integration in time (or, equivalently, in azimuth \j/) until the following 
convergence criteria are met: 

s[2(n + l)n/Q] = s[2nJt/Q] ± £i 
i[2(n + 1 )jc/Q] = j[2«7t/£2] ± 82 


where 


n = number of rotor rotations 
Ej, e 2 = convergence tolerances (m, m/s) 
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R 

M= J my 2 d ^ 

o 

R 

Mr = J mrydq 

o 

R 

M b = J e-ntnydl, 

o 

R 

M g = J mydi, 

o 

R 

Kb= Je/ C (y"M 
o 

R 

Fa = J 7h(^)(y ') 2 d\ 

o 

R 

K q = I m (R)y'(R)y(L) - J I m (y’) 2 DE 

0 

r 

K c = J W(y 'fd% 

o 

R 

F a = J dAydl 

0 

R 

Tail) = J mrd \ 

% 

RT 

T c d ) = J ™yd\ 

5 

R 

Tgd) = J md\) 


( 11 - 32 ) 


After a trim solution has been found or after each revolution of a yaw-motion solution, the 
flap displacement, slope, and velocity can be computed from the following equations: 


V = s(t)y(£) 

(ll-33a) 

v' = s(t)y'd) 

(1 l-33b) 

V = s(t)y(E,) 

(ll-33c) 
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Rotor Blade Loads 

Blade loads distributed along the ^-axis are calculated by the force-integration method, 
as follows: From Equations (11-27) and (ll-32c), the spanwise tension force load is 

m)= n 2 T a ({,)+2ClsQ p sT c (?,)-g C yT g a) (11-33) Q3 

Integrating Equation (ll-4b) gives the blade flapwise shear force load as 

R 

V C (5) = J dA^- (O 2 p 0 + 2Q^cv|/ + ^V|/)c0 p 7h(^) 

S 

+ (O 2 s0 2 s - s)T c (q) + n 2 sQ p cQ p M e Cq) (11-34) 

+ g(~X cQ p + sQ p s\\i + Poc0 p cv|/)r g (^) 

Integrating Equation (ll-4a) gives the blade chordwise shear force load as 

R 

Vt,(£) = J dA T1 ^ + (Q 2 p o + 2Q<i)c\|/ + ^t|/)^0 p 7n(^) 

I 

+ Q 2 sQ p cQ p sT c (i,) + Q 2 c0 2 M e (£) (11-35) 

+ g(X*% + cQpsy - PosQpC\\f)T g (f,) 

Integrating Equation (ll-5c) gives the blade pitch moment load as 

R R 

M^) = -jv^)V(&dZ, + J qU^dl 

\ 5 (11-36) 

- C1 2 sQ p cQ p AI ( ^) + gs\\tsQpM e (fy 


Integrating Equation (ll-5b) gives the blade edgewise moment load as 


R R 

M&) = jv^)d^-n 2 je^)rdl 

5 5 (11-37) 

R 

-Q 2 cQpsQ p s \l^l)dl + gcyMfff) 
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Finally, integrating Equation (1 l-5a) gives the blade flatwise moment load as 

R R 

M n (£) = - \vtf£)d\+ 

Z, k 

R 

k 

In the above equations, 

R 

M&) = { e^)m{l)dl 

Z, 

R 

AZ(^) = J - I m md% 

% 


(11-38) 


(11 -39a) 
(11 -39b) 


Typical Computer Model for Blade Load Calculations 

The calculation of wind turbine blade loads for even the simplest of cases, such as the 
one described here ( i.e ., one degree of freedom and one blade mode), requires a computer 
model. A variety of special-purpose models for HAWT loads have been developed in the 
past, dating back to the 1970s when only mainframe computers had sufficient capability for 
the task [Spera 1977]. Today, HAWT load computer models for personal computers are 
available, such as the FLAP Code [Wright et al. 1988], The following discussion is based 
on experience with the FLAP Code, which is a current analysis tool in the public domain. 

Requirements 

A computerized solution of the equation of motion and computation of displacements 
and loads requires a sophisticated, interactive model that is flexible, well-documented, and 
easily modified. The model must be capable of performing a variety of tasks, including 

- input of data and output of results; 

- matrix inversion; 

- time-domain analysis of differential equations; 

- computation of spatially-dependent blade properties and aerodynamic parameters. 

The computational requirements of the equation of motion and other quantities associated 
with it determine which portions of the model can forego efficiency in favor of flexibility 
and readability, and which portions need to be as efficient as possible. This can be 
accommodated in a computer model composed of two modules, as illustrated by the flow 
chart in Figure 1 1-4. 

Module 1 is a data pre-processor, in which a file of raw blade and turbine property is 
processed to produce a data file that can be used to solve the equation of motion. 
Flexibility and readability are more important than efficiency in this module. Module 1 also 
computes all the coefficient matrices since they are independent of most of the non- 
structural variables, such as wind speed and tower shadow. Module 2 performs the actual 
model run including solution of the equation of motion, computation of the loads and output 
of the results. Here, efficiency is a primary requirement. 
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Module 1 


Module 2 




Figure 11-4. Typical flow chart for a computer model used to calculate HAWT blade 
loads. [Thresher et al. 1986] 


Computational Procedure 

The procedure followed in the FLAP Code for solving Equation (11-31), the 
fundamental equation of motion for a HAWT blade, can serve as an illustration of an 
approach that has been used successfully by the creators of such computer models. The 
general approach is as follows: 

1. All coefficient matrices are multiplied by the inverse of the mass matrix. Thus the 
mass matrix associated with blade tip accelerations multiplied by its inverse gives the 
identity matrix. This results in a set of equations with the blade tip accelerations on 
one side of the equality and the computed generalized forces on the other. Multiplica- 
tion by the inverse mass matrix is only done once, at the beginning of Module 2. 
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2. In this way, the accelerations associated with each flap coordinate function can be 
evaluated numerically be substituting the current values for the flap velocities and 
displacements into the force side of the equation. 

3. The blade tip displacements are computed by solving the second-order differential 
equation relating acceleration, velocity, and displacement. The solution is performed 
via the modified Euler predictor-corrector method, which uses the current blade tip 
accelerations and the previous values of the displacements and velocities. 

4. The blade loads are computed only at the completion of a trim solution. 

5. The solution during time-dependent prescribed yawing motions is run at the completion 
of the trim solution. Loads are computed at the completion of each revolution during 
the yawing solution. Yawing motion is prescribed according to the simple equation 

<|>(0 = <(>o+ itJflSinCGfyt) (11-40) 

6. A set of nine blade motion and blade load parameters is computed for each specified 
spanwise location on the blade axis (typically 1 1 stations) at each azimuthal position 
of the blade during a revolution (typically 36 azimuths). These parameters are the 

- flap displacement, v; 

- flapwise slope, v'; 

- flap velocity, dv/dt, 

- edgewise shear load, 

- flapwise shear load, 

- spanwise tension load, T; 

- flapwise moment load, 

- edgewise moment load, M ? ; 

- pitching moment load, AL. 

Sample Blade Load Analysis 

The turbine modeled in this example is the Grumman prototype HAWT (see Table 3-1), 
which has a three-bladed, rigid-hub rotor 10.1-m in diameter that runs downwind of a shell 
tower [Adler et al. 1980]. The lift and drag models are simplifications of the airfoil data 
presented by Sweeney and Nixon [1978]. Other input data are as follows: 


U H 

= 9.1 m/s 

m 

= 0.143 

1 0 

= 0.25 

t 

= 0 

Vo 

= 15 deg 

♦o 

= 10 deg 

Q 

= 7.54 rad/s 

c 

= 0.46 m 

Po 

= 3.5 deg 

e, 

= 0 deg 


The natural frequency of the blade’s first ( i.e . lowest frequency) flapwise bending mode is 
approximately 3.95 Hz, or about 3.3 times the rotor rotational frequency. 

Flap bending moments in the operating turbine were measured at a spanwise station 1.8 
m from the rotor axis (20% of span) and averaged over 20 revolutions. The test data were 
azimuth-averaged in this manner to remove random wind effects that are not modeled in 
this analysis. After averaging, the mean value of the measured flap bending moment was 
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subtracted, producing the variable flap moment history labeled “Test Data” in Figure 11-5. 
The calculated flap moment load at 20% of span, again with the average removed, is also 
plotted in Figure 11-5, permitting a comparison of the blade’s flapping response to gravity, 
tower shadow, and average wind shear. 

The test and calculated wave forms are very similar in Figure 11-5, indicating a 
validated computer model for this turbine. However, another HAWT would require a 
similar test/calculation comparison for validation of its structural, aerodynamic, and wind 
input data. These data are an integral part of the blade load analysis model, and as critical 
to validation as the basic equations of motion. 



Figure 11-5. Sample comparison between calculated and measured flapwise moment 
loads. The test turbine is the 10.1-m diameter Grumman prototype HAWT. [Thresher et 
al. 1986] 


Current and Future Developments in Structural-Dynamic Modeling 

The demands placed on structural dynamics models are becoming more severe as 
turbine designs evolve. Models are being asked to account for more-complex phenomena, 
including 

— non-steady air flow; 

— blade-to-blade interactions; 

— yaw and pitch motions of the nacelle and tower; 

— blade pitch control; 

— teetering motion; 

-- variable-speed; 

— startup transient; 

— shutdown and braking transients. 


ASM_Wind Turbine_Ch1 1 .indd 


632 


MTC 


02/16/2009 10:06AM 



Structural Dynamic Behavior of Wind Turbines 633 

The desire to account for all such effects within one or two comprehensive models is 
driving engineers to develop computer codes of ever-increasing complexity. 

Improved structural dynamic models of wind turbines appear to be developing along 
two different paths. The first approach involves special purpose codes, each written to 
analyze the structural dynamic behavior of a specific configuration of wind turbine. Such 
codes generally have 12 to 15 degrees of freedom (DOF) tailored to the significant motions 
of the “target” structural system. The second approach seeks to adapt general purpose 
structural analysis codes to simulate wind turbine dynamics. These general purpose codes 
are powerful enough to handle almost any turbine configuration, but they are complex and 
expensive. Moreover, the structural dynamicist must have a significant level of expertise 
in their use, and each new configuration requires considerable time for modeling and 
verification efforts. 

Typical Special-Purpose Model 

The structural analysis program PHATAS-II, developed at the Netherlands Energy 
Research Foundation (ECN), is an example of a current special-purpose code [Lindenburg 
and Snel 1993]. The PHATAS computer model can be used to calculate dynamic displace- 
ments and loads for the main structural subsystems of a HAWT. Wind inputs to this model 
include 

— steady wind, including vertical and horizontal gradients; 

— wind direction as a function of elevation; 

— rotationally-sampled wind turbulence; 

— turbulent wind in time-series form. 

A number of structural degrees of freedom are available within the code, such as 

— flapwise and chordwise elastic blade deformations; 

— blade pitch motion; 

— blade flapping and teetering hinges; 

— torsional deflection of the power train; 

— variable rotor speed; 

— pitch and yaw motions of the nacelle; 

— bending and torsional deflections of the tower. 

Typical General-Purpose Model 

The ADAMS code, developed at Mechanical Dynamics, Inc., is an example of a 
general-purpose, multi-body structural analysis program that has been adapted for wind 
turbine applications [Malcolm and Wright 1994], This code has virtually unlimited degrees 
of freedom and places no limits on the type or magnitude of displacements. While the 
ADAMS code can be run on a high-end personal computer (PC), it is exceptionally slow 
when modeling systems with many degrees of freedom or when performing simulations 
involving wind turbulence inputs. For example, it can take several days of PC run time to 
simulate 10 minutes of real time operation of a HAWT subjected to turbulent winds. 

The strengths of this type of computer model are that almost any structural system 
(wind turbine or otherwise) can be modeled with it once the required skills have been 
learned by the user, and a mixture of complex external loadings and operating conditions 
can be applied. As PC performance increases, this type of general-purpose computer model 
is expected to used more and more for the structural analysis of modem wind turbines. 
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Natural Vibration Modes in a Wind Turbine 


Dynamic effects can be substantial in a wind turbine system because of the periodic 
nature of its aerodynamic loading. As a result, the design process must depend on accurate 
predictive tools and/or measurement techniques for determining the dynamic response 
characteristics of the rotating structure. An elastic structure responds to periodic forcing 
functions by vibrating in one or more discrete geometric patterns called mode shapes, each 
of which has a companion rate of vibration called a modal frequency. Mode shapes and 
frequencies are determined primarily by the distributions of mass and stiffness throughout 
the structure and by its boundary conditions. Rotation can alter the natural frequencies of 
certain mode shapes, when centrifugal and Coriolis forces change stiffnesses. Because 
mode shapes and frequencies are relatively independent of applied loadings, they are often 
referred to as natural or free-vibration modes. Both analytical and experimental procedures 
for determining mode shapes and frequencies are referred to as modal analysis. 

Figure 11-6 shows the results of a typical modal analysis of a HAWT [Sullivan 1981], 
displayed in a frequency plot called a Campbell diagram. Natural frequencies are plotted 
vs. rotational speed, with one generally-horizontal line per mode. Rays from the origin are 
plots of integer multiples of the rotational frequency. In a rotating structure, excitation 
loadings normally occur at these integer-multiple frequencies, often abbreviated IP, 2 P, etc. 


Mode: 



Rotor Speed (rpm) 


Figure 11-6. Typical Campbell diagram displaying the results of a HAWT modal 
analysis. Radial lines indicate the frequencies of potential excitation forces [Sullivan 1981] 
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Defining modal characteristics is of paramount importance in that, with this 
information, potentially harmful resonant behavior can be avoided. Resonance is a 
phenomenon which occurs in a structure when the frequency of periodic loading or 
excitation equals or nearly equals one of the modal frequencies of the structure. Thus, in 
Figure 11-6, every intersection of a radial line and a modal-frequency line is a potential 
resonance. Resonance if often characterized by a large increase in displacements and 
internal loads. The relative size of this increase is measured by a dynamic amplification 
factor, which is the ratio of the amplitude of the dynamic displacement to the static 
displacement under the same forcing function. 

Not all resonances cause large amplification factors, because rotation can sometimes 
eliminate the matching of a forcing function with a natural mode, even when both have the 
same frequency in a non-rotating frame of reference. A general discussion of amplified and 
non-amplified resonances is beyond the scope of this chapter, but background information 
on resonances significant in HAWT systems is given in [Sullivan 1981]. 

In rotating structures, the acquisition of dynamic response characteristics is complicated 
by their dependence on the structure’s rotational speed. Therefore, special methods have 
to be employed for both the predictive and experimental modal analysis. Predictive 
techniques for wind turbines are usually based on the finite-element method, in order to take 
advantage of its superior versatility. As the turbine is modeled in a frame of reference 
which rotates at the turbine operating speed (assumed constant), rotating coordinate system 
effects must be taken into account. 

The finite-element method has been used extensively to predict the modal characteris- 
tics of rotating structures. When this method is applied to HAWT blades and other 
structures which are comprised of rotating beams reminiscent of fan or helicopter blades, 
some of the rotating coordinate system effects are neglected. This is usually justified when 
motions perpendicular to the axis of rotation {i.e., along the axis of a HAWT blade) are 
relatively small. This is not the case for a VAWT, however, and these effects cannot be 
omitted. Here we will use a method similar to the one described by Patel and Seltzer 
[1971] for analyzing the modes of the spinning “Skylab” satellite, in which all rotating 
coordinate system effects are included. 

Other analytical techniques which have been used to predict the dynamic characteristics 
of rotating structures include various Galerkin procedures [Klahs and Ginsberg 1979, 
Warmbrodt and Friedmann 1980], the modal method, and the component mode method in 
combination with various eigensolution schemes [Vollan and Zwaan 1977, Ottens and 
Zwaan 1978a and 1978b]. A major experimental obstacle to obtaining the natural 
frequencies and mode shapes of a rotating structure is that mechanisms must be devised for 
exciting the structure while it is rotating. Sliprings or telemetry must also be employed for 
data transmission from instruments such as accelerometers on the structure. 


Finite Element Modal Analysis of a VAWT 

When performing a modal analysis of a rotating structure, the motion is most 
conveniently measured relative to a coordinate system that is rotating with the structure. 
This permits the displacements to be small by eliminating the large rigid-body rotation. The 
coordinate system employed in this analysis rotates at the turbine operating speed, Q, with 
the origin fixed in space at the lower bearing of the rotor. The angular velocity vector is 
directed vertically upward along the positive z axis. The finite-element matrices are 
developed on the basis of the motions that remain small within this rotating frame. For an 
accurate representation of the dynamics of a VAWT relative to a coordinate system rotating 
at a constant speed, the following effects must be included: 
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— tension stiffening', 

— centrifugal forces', 

— Coriolis accelerations. 

The tension stiffening that affects primarily the blades is induced by steady centrifugal and 
gravity forces. The centrifugal and Coriolis terms are a direct result of using a rotating 
coordinate system. 

Equations of Motion 

The total acceleration at a point in the rotating structure, with respect to fixed 
coordinates, can be represented in terms of the acceleration in the rotating coordinate system 
by an equation similar to Equation (11-18), as follows: 

a t = u + IClxu + Qx[Qx(r + a)] (11-41) 


where 

a t = total acceleration vector, excluding gravity (m/s 2 ) 
r = original position vector of an element; origin at lower bearing (m) 
u = displacement vector; observed in the rotating coordinate system (m) 

Q = constant angular velocity vector of the rotating coordinate system (rad/s) 

Using this expression for the acceleration in the equation of motion for a structure, the usual 
damping and stiffness matrices are altered from those of a static structure. The resulting 
differential equations are represented by 


Mu + Cu- Su + Ku = F c + F g 


(11-42) 


where 

M 

u 

C 

s 

K 

F 

C 

F 


normal (unaltered) mass matrix (kg) 

displacement observed in the rotating coordinate system; dot signifies 

derivative with respect to time (m) 

additional Coriolis matrix (N-s/m) 

additional centrifugal softening matrix (N/m) 

normal (unaltered) stiffness matrix (N/m) 

additional static load vector representing steady centrifugal force (N) 
gravitational forces; steady in a VAWT (N) 


The Coriolis matrix, C, is skew-symmetric and results from the second term on the right 
side of Equation (11-41). Here we have ignored any structural (internal) damping because 
most VAWTs are lightly damped. Therefore, the C matrix is totally the result of Coriolis 
effects. The centrifugal softening matrix, S, comes from the variable portion (i.e., 
dependent on k) of the last term on the right in Equation (11-41). F c also comes from the 
last term in Equation (11-41). To obtain the mode shapes and frequencies of the turbine 
as observed in the rotating system, Equation (11-41) is reduced to the following form: 


Mu + Cu + (K+ K g - S) u = 0 


(11-43) 


where K ( . = geometric stiffness matrix resulting from steady centrifugal and 
gravitational loadings (N/m) 


ASM_Wind Turbine_Ch1 1 .indd 


636 


MTC 


02/16/2009 10:06AM 



Structural Dynamic Behavior of Wind Turbines 


637 


Thus, the solutions correspond to small motions about a prestressed state. Aeroelastic 
effects ( i.e ., aerodynamic loading caused by structural motions) have not been included in 
this analysis. While strongly influencing the aerodynamic stability of a VAWT [e.g., 
Popelka 1982], aeroelasticity seems to have little effect on the natural frequencies and 
modes of the turbine. 

Eigenvalue Solution 

Obtaining mode shapes and frequencies from Equation (11-43) involves solving what 
is known as a characteristic-value problem, which has the standard form 


011*1 + 012*2 + — + 0 u x n = CO*! 

<* 21*1 + 022*2 + — + 02 n x n = “*2 


(11-44) 


0 « 1*1 + 0 « 2*2 + ■■■ + a nn x n = “*, 


where, for our modal analysis problem 

n = number of degrees of freedom; number of finite elements times three 
displacement components per element 
X ,„x n = eigenvector or characteristic vector, defines one mode shape 

co = eigenvalue or characteristic value-, defines one modal frequency (rad/s) 

Each non-trivial solution of this characteristic-value problem has an eigenvector and a 
corresponding eigenvalue, and the set of eigenvector/eigenvalue pairs comprises the results 
of a theoretical modal analysis of our structural system. Modes are usually presented in 
ascending order of frequency, from the lowest or fundamental mode to the highest fre- 
quency of interest. The latter seldom needs to exceed ten times the maximum rotational 
speed of the turbine (10P). 

To find the eigenvalues of the system of equations represented by Equation (11-43), 
it is more convenient to use state vectors instead of the displacements u. Hence, we will 
introduce the state vector 


u 


V = 


u 


and the matrices 


M | 0 

--I-- 

. 0 | X'. 
C I I? ' 

--I-- 

-K 1 | 0 


(11 -45a) 


(11 -45b) 
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in which 

K'= K+K G - S 

Equation (1 1-43) can now be written as 

Air + B v = 0 (11-46) 

Here A is real symmetric, and B is real skew -symmetric. Seeking a solution for Equation 
(1 1-46) in the form 

v(t) = $e im (11 -47a) 

where i is the square root of -1, the resulting eigenvalue problem is 

\B + iu)A](j) = 0 (11 -47b) 

We now define an eigenvector x as 


x = A m i j> (11 -47 c) 

assuming that A is now positive definite. We can now transform Equation (11 -47b) into 
the standard eigenvalue form of Equation (11-44), with all terms on the left side, as follows: 

[G-oo/]x=0 (11 -48a) 

G = iA ll2 BA 112 (11 -48b) 

Because of the skew-symmetry of the matrix B, the matrix G is Hermetian. Consequently, 
it can be shown that the eigenvectors are, in general, complex; but the eigenvalues are real 
[Franklin 1968]. 2 If structural damping is included, the system loses it Hermetian character 
and the eigenvalues as well as the eigenvectors are complex. Again, in the present analysis, 
structural damping has been ignored because VAWTs are lightly damped. The Hermetian 
character of the eigenvalue system has important ramifications to the rotating modal test. 
Even if the modes of the non-rotating structure are perfectly real, the modes of the rotating 
structure will be complex. 


2 If A is not positive-definite, then G is no longer Hermetian and the eigenvalues are 
not necessarily real. This condition can lead to dynamic instability. 
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NASTRAN Finite Element Modeling 

The general-purpose structural analysis code NASTRAN is the software commonly used 
for modeling a wind turbine in finite elements to determine its natural vibration modes. 
Before presenting an example of the modeling, a short description of the NASTRAN calcu- 
lational procedure is in order, as it is applied to this specific problem. Two of the available 
NASTRAN rigid format options are used in the modeling effort. First, a static analysis 
with geometric stiffening effects is performed on the model under the action of centrifugal, 
gravitational, and boundary forces. The resulting modified stiffness matrix ( K + K a in Eq. 
(11-43)) is then retained via a modest amount of DMAP programming for use in the 
subsequent complex eigenvalue analysis (the second rigid format). Along with necessary 
structural data, the resulting data deck includes so-called “direct matrix input at grid points” 
( DMIG ) cards. It is through these that the Coriolis and centrifugal softening matrices (C 
and S in Eq. (1 1-43)) are input to the model. 

Because the required number of DMIG cards tends to be large, a pre-processor is 
generally used to make the inclusion of rotating coordinate system effects in NASTRAN 
relatively transparent to the user. This pre-processing program 

— reads the NASTRAN structural data deck and extracts needed information from 
GRID, CBAR, MAT1, CONM2, and RFORCE cards; 

— constructs the necessary matrices and casts them in the DMIG card format; 

— inserts the DMIG data cards into the NASTRAN data deck. 

NASTRAN users are provided with a choice of several complex eigensolvers. One 
employs the upper Hessenberg method and is very efficient. A method that has been used 
successfully for the modal analysis of wind turbines is to first model the structure with the 
necessary number of degrees of freedom, then use a standard Guyan reduction procedure 
(also a NASTRAN option) to reduce the degrees of freedom to a manageable number 
(approximately 100) and then apply the upper-Hessenberg method. 
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Sample Theoretical Modal Analysis 

In this sample modal analysis, taken from [Came et al. 1982], modal shapes and fre- 
quencies are determined for the Sandia/DOE 2-m VAWT shown in Figure 11-7. The 
configuration of the finite-element model of this VAWT is shown in Figure 11-8. The 
model is composed of two blades, a central column, and elastic supports. The curved 
blades are modeled with 20 beam elements each. The truss column, which adds a fair 
degree of complexity to the model, is represented by 150 beam elements. The masses of 
the hardware associated with the upper bearing and its support-cable connections are repre- 
sented by a single concentrated mass placed at the top of the truss model. The support 
cables are modeled with horizontal, linear springs and a vertical downward force on the 
truss, as shown in Figure 11-8. 

The masses of the test instrumentation installed on the rotor (which are not negligible 
compared to the mass of the blades) are included by adding appropriate concentrated 
masses. The base structure is modeled as a concentrated mass in combination with torsional 
and linear springs. The horizontal springs at the upper and lower bearings have equal 
stiffnesses in two orthogonal directions, producing the desired isotropic boundary conditions 
typical of VAWTs. The stiffness and mass properties of the support system must usually 
be estimated, at least to some degree. The initial estimates for this model are given in 
parenthesis in the figure. The final properties of the supports were obtained by “tuning” the 
model to match frequencies measured when the VAWT was not rotating, as will be 
discussed in the next section. 
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Figure 11-7. Sandia/DOE 2-m test VAWT with typical modal-survey instrumentation 
and excitation devices. ( Courtesy ofSandia National Laboratories ) 
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1 

2 

3 



26 

27 

Non-Rotating System Modes 28 

29 

Table 11-1 lists the frequencies of the first ten modes of vibration on the non-rotating 30 

C i.e ., parked) 2 -m VAWT, with its rotor brake on. The first column gives the mode 31 

number, usually assigned in increasing order of frequency. Mode identification terminology 32 

is given in the second column. The corresponding mode shapes, in three views, are 33 

sketched in [Came et al. 1982 ], The third column lists the frequencies measured with the 34 

rotor parked, and the fourth column gives the frequencies calculated using initial estimates 35 

of the support mass and stiffness properties. Again, these estimates are the data in 35 

parenthesis in Figure 11 - 8 . In the fifth column of Table 11-1 are the “tuned” calculations 37 

of the modal frequencies for the non-operating condition. The last column in the table 3g 

gives the relative changes in frequency resulting from tuning the model. 39 

40 

Tuning the Finite-Element Model 41 

42 

The differences between the initial and tuned frequencies, ranging from - 20 % to +8%, 43 

indicate the scale of the errors that can be present when (1) the support system and the 44 

resulting theoretical boundary conditions are complicated, and (2) some of the auxiliary 45 

components in the structural system are not modeled in detail. For practical reasons, the 46 

latter is usually the case. As a result, a finite-element model must be calibrated against a 47 

basic set of measured responses, on the same or a similar structure, to establish confidence 4g 

in predictions made with it of natural frequencies under operating conditions. Modal 49 

analysis is clearly recognized to be a combination of theoretical and experimental methods. 50 

51 
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Table 11-1. 


Modal Frequencies of the Non-Rotating 2-m VAWT [Came et al. 1982] 


Mode 

Number 

Mode name 

Test 

(Hz) 

Initial 

model 

(Hz) 

Tuned 

model 

(Hz) 

Tuning 

change 

(Hz) 

1 

1st Anti-Symmetric Flatwise 

12.3 

12.5 

12.3 

-1.6% 

2 

1st Symmetric Flatwise 

12.5 

12.6 

12.4 

-1.6% 

3 

1st Rotor Out-of-Plane 

15.3 

17.1 

15.2 

-11.1% 

4 

1st Rotor In- Plane 

15.8 

17.2 

15.9 

-7.6% 

5 

Dumbbell 

24.4 

22.6 

24.4 

+8.0% 

6 

2nd Rotor Out-of-Plane 

26.2 

30.5 

26.2 

-14.1% 

7 

2nd Rotor In-Plane 

28.3 

30.6 

28.0 

-8.5% 

8 

2nd Symmetric Flatwise 

29.7 

30.9 

30.6 

-1.0% 

9 

2nd Anti-Symmetric Flatwise 

31.5 

39.7 

31.7 

-20.2% 

10 

3rd Rotor Out-of-Plane 

36.5 

42.3 

36.5 

-13.7% 


Calculated Rotating System Modes 

The Campbell diagram in Figure 11-9 presents the modal frequencies of the 2-m 
VAWT at rotor speeds from 0 to 600 rpm. The curves in this figure are the calculated 
frequencies obtained using the tuned NASTRAN model, with mode numbers corresponding 
to those in Table 1 1-1. The variations of the frequencies with rotor speed are quite complex. 
While most increase monotonically with speed, some decrease monotonically. Others 
increase and then decrease and vice versa. These variations are similar to those associated 
with the classical whirling-shaft problem. In addition to complicated frequency variations, 
the mode shapes also change in character with rotor speed. Specifically, mode shapes 
which are completely uncoupled when the turbine is parked become coupled during rotation. 
This coupling is discussed in more detail in [Came et al. 1982], In mathematical terms, the 
modes change from real to complex as the turbine rotates. 


Sample Experimental Modal Analysis 

An overview of the general techniques for modal testing of wind turbines will now be 
presented, again using the Sandia/DOE 2-m VAWT as an example. Portable data acquisi- 
tion tools of the test engineer are based on the fast Fourier transform. The FFT technique, 
which involves exciting the structure with a force having a linear spectrum containing the 
frequency band of interest, is generally faster and more versatile than the classical swept- 
sine technique. The applied forces and responses are measured in the time domain and 
transformed to the frequency domain using the FFT. The frequency response functions are 
then computed from the cross-spectral and auto-spectral densities of the applied force and 
the responses. 

Typically, several measurements are averaged to reduce the effects of uncorrelated 
noise. A more-complete description of FFT modal testing is contained in [Klosterman and 
Zimmerman 1975]. The greater versatility of this technique is a result of its more-relaxed 
requirements on the excitation force. For example, the FFT technique is equally applicable 
to shaker-driven excitations, instrumented-hammer impacts, and excitations by wind forces. 
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Figure 11-9. Campbell diagram of the Sandia/DOE 2-m VAWT, presenting both theo- 
retical and experimental modal frequencies vs. rotor speed. [Came et al., 1982] 


Instrumentation 

The general layout of instrumentation on the test VAWT is illustrated in Figure 1 1 -7. 
For non-rotating tests, the turbine was instrumented with 18 piezoelectric accelerometers 
on the blades, 10 on the central column, 20 on the supporting pedestal, and 6 on the turbine 
shaft assembly. All input and response signals were low-passed at 100 Hz, amplified as 
necessary with signal-conditioning amplifiers, and recorded on a 14-channel tape recorder. 
In the rotating modal test the high, steady centrifugal accelerations prevented the use of 
accelerometers except on the central column. Strain gages are normally installed on newly- 
designed wind turbines undergoing development tests, and these gages are commonly used 
as response transducers for modal tests as well. 
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For this rotating test, strain gages were placed on the blades and central column, using a 
double-active gage in the bridge so that bending about a single axis would be measured. 
Gages were placed near the blade roots on leading and trailing edges (to measure chordwise 
bending), on the inside and outside surfaces at maximum blade thickness (to measure flat- 
wise bending), and on the central column (to measure in-plane and out-of-plane bending). 
On-board amplifiers boosted the strain signal to 2,000 psi/V for all gages, in order to give an 
adequate level for passing through the sliprings and the long wires that are typically required 
for modal tests on operating wind turbines. Because of expected difficulties in understanding 
the modes of the rotating central column, two piezoresistive accelerometers were attached to 
the upper column to record in- and out-of-plane motions. 

Excitation 

The scope of methods used to excite vibrations in wind turbines for modal analysis is 
broad. These methods fall into the four general categories of human, impact, step-relaxation, 
and natural-wind. As mentioned previously, the FFT-testing procedure allows flexibility in 
the excitation method, which is of significant benefit. 

Human 

Human excitation is the simplest approach to implement, but has the most limitations. 
The large displacement, low frequency modes of large wind turbines allow a test engineer to 
manually input a quasi-sinusoidal force into the structure to excite a single mode of the sys- 
tem. The frequency and damping parameters of the mode can be determined as the systems 
rings down after the excitation is removed. This approach was used to test over 20 commer- 
cial 17-m VAWTs in as little as 1 V* hours each. This study also provided valuable insight 
into the structural-dynamic variability of duplicate units in the field [Lauffer, 1986]. 

Impact 

Impact excitation methods range from using an instrumented 3-lb hammer to strike the 
2-m VAWT at various locations on its blades and central tower to swinging a 1 ,000-lb ram 
from a crane to strike the tip of the NASA/DOE 2.5-MW Mod-2 HAWT [Boeing 1982], In 
both cases, the impact force must be of sufficient amplitude and duration to excite the modes 
and frequencies of interest. This is accomplished by prior calibration tests. 

Impact testing was also one of the techniques used to test the Sandia 34-m VAWT Test 
Bed in various stages of assembly. This allowed model correlation of substructure models 
before final model assembly. The substructures tested included two blade sections and a 
blade assembly. Step-relaxation and natural wind excitation were also used for other sub- 
structures, as appropriate [Came etal., 1989]. 

Impact testing was also used to study the use of a Laser Doppler Velocimeter (LDV) 
for non-contact measurements of a HAWT in the field. This study used impact excitation to 
compare the LDV to traditional accelerometers as well as to study damage detection tech- 
niques. A comparison to natural wind excitation was also performed. The test object in this 
study, a 15 m diameter HAWT in its parked configuration [Rumsey et al., 1997]. 
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Step -Relaxation 

Step-relaxation relies of the quick release of strain energy in a pretensioned cable for 
the excitation of the structure. A small device applying a force of approximately 50 lb was 
used successfully on the 2-m VAWT [Came et al. 1982], A heated wire was used to bum 
through a nylon cord to affect the quick release of the excitation force. For a modal test of 
the 64-m diameter Eole (Fig. 3-38) two separate force-application points were used, one on 
the central column and one on a blade [Came et al. 1988]. Forces were applied through high- 
strength steel cables loaded with a diesel-powered winch and restrained by a nylon strap. 
Loads of 135 kN (30,000 lb) and 45 kN (10,000 lb) were applied to the column and blade, 
respectively. 

One difficulty in using step relaxation is related to the FFT algorithm. A quick release 
basically applies a force which is a Heaviside function, for which there is no Fourier trans- 
form. However, if the Heaviside function is passed through a high-pass filter, it can be con- 
verted to a well-behaved function that is easily transformable with the FFT. 

Natural Wind 

Using the natural wind to excite the turbine has certain advantages, as discussed in 
[Came et al. 1988]. Two of the most significant are the greatly reduced cost and complexity 
of the test and the ability to test during windy conditions. Wind-excitation methods were 
developed during test of the Eole VAWT, and these produced modal frequencies in excellent 
agreement with those obtained by step relaxation. However, damping information is not as 
readily available from the power spectra obtained using wind excitation. Development of 
the Natural Excitation Technique (NExT) rectified this situation and will be discussed below 
[James etal., 1995]. 

Measured Rotating System Modes of the 2-m VAWT 

The results of the modal testing performed on the 2-m VAWT using step-relaxation 
excitation are shown by the data points in Figure 11-9. The correlation here is typical of 
wind turbine modal analysis. The average absolute deviation of the experimental frequen- 
cies from the theoretical (tuned-model) frequencies is 0.5 Hz, or 2.2%. Deviations generally 
increase with increasing rotor speed, which is a characteristic of finite-element models that 
do not have a sufficient number of degrees of freedom to accurately represent complex mode 
shapes. 


ASM_Wind T urbine_Ch1 1 .indd 


645 


MTC 


02/16/2009 10:06AM 



646 


Chapter 11 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 
11 
12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 


Development of Natural Wind Excitation Testing 

The structural dynamics of rotating VAWT machines are critical design conditions that 
require experimental techniques that can match the fidelity of the advanced analyses described 
previously. The previous example using the 2-m VAWT showed the power of performing 
step-relaxation testing during rotation. However, this approach is technically intensive and 
would not be generally applicable for many operational scenarios. As mentioned previously, 
the use of natural wind excitation is a more generally applicable approach and allows test- 
ing in the actual operating environment. This section details the development of the natural 
excitation technology in the form of the natural excitation technique or NExT [James et al., 
1995]. A more complete background on the historical development of the natural wind exci- 
tation testing for wind turbines is also available [Came and James 2008]. 

Theoretical Development of NExT 

A critical step in the development of NExT was to find a relationship between modal pa- 
rameters of structural interest (such as frequency and damping) and a function that could be 
produced from measured response data. For NExT, the cross-correlation function between 
two response-only measurements was used. The derivation of the relationship of modal 
parameters and the cross-correlation function begins by assuming the standard matrix equa- 
tions of motion: 

[M][3c(m + [C][i(0] + [*]{*«) = (1 1-49) 

where 

[M\ = mass matrix (g) 

[C] = damping matrix(N-s/m) 

[AT] = stiffness matrix (N/m] 

[f] = vector of random forcing functions (N) 

[jc] = vector of random displacements (m) 
t = elapsed time (s) 

x, x = first and second derivatives with respect to time (m/s, m/s 2 ) 

Equation (11-49) can be expressed in modal coordinates using a standard modal transfor- 
mation after performing a standard matrix diagonalization (assuming proportional damping). 
A solution to the resulting scalar modal equations can be obtained by means of a Duhamel 
integral, assuming a general forcing function, {/}, with zero initial conditions (i.e. zero initial 
displacement and zero initial velocity). The solution equations can be converted back into 
physical coordinates and specialized for a single input force and a single output displacement 
using appropriate mode shape matrix entries. The following equations result: 

Xikit) = ^VirWkr • f f fk(t)g r (t - l)dx (11-50) 

T= 1 

g r (t) = 0, for t < 0 

g r (t) = eX P sin f° r t^0; 

where 

®d = ®n (1 - tf 2 y n is the damped modal frequency (rad/s) 

G> = rth modal frequency (rad/s) 


ASM_Wind Turbine_Ch1 1 .indd 


MTC 


02/16/2009 10:06AM 



Structural Dynamic Behavior of Wind Turbines 


647 


£ r = rth modal damping ratio 
m r = rth modal mass (g) 
n = number of modes; 

\|r. = ith component of mode shape r 

The next step of the theoretical development is to form the cross-correlation function of 
two responses (x k and x jk ) to a white noise input at a particular input point k. The cross- 
correlation function 9i. jk (7’) can be defined as the expected value of the product of two re- 
sponses evaluated at a time separation of T: 

^ijk(T) = E[x ik (t + T) xjkit )] (11-51) 

where E is the expectation operator. 

Substituting Equation (1 1-50) into (11-51), recognizing that f f t) is the only random vari- 
able, allows the expectation operator to depend only on the forcing function. Using the defi- 
nition of the autocorrelation Junction, and assuming for simplicity that the forcing function 
is white noise, the expectation operation collapses to a scalar times a Dirac delta function. 
The Dirac delta function then collapses one of the Duhamel integrations embedded in the 
cross-correlation function. The resulting equation can be simplified via a change of vari- 
able of integration, setting X = M. Using the definition of g from Equation (11-50) and the 
trigonometric identity for the sine of a sum results in all the terms involving T being separated 
from those involving X. This separation allows terms that depend on T to be factored out of 
the remaining integral and out of one of the modal summations, with the following results: 

„ Aijk exp (- C < t) cos t) 

Xtjk(T) = t (H-52) 

r=1 + B. r . k exp T J sin 

where A' jk and B'. k are independent of T, are functions of only the modal parameters, and 
completely contain the remaining modal summation, as shown in the following equations: 



£ «£\|r»-\| tkrVjsVks 
* =1 m r i ifj 





sin 




(11-53) 


Equation (11-52) is the key result of this derivation. Examining Equation (11-52), the 
cross-correlation function is seen to be a sum of decaying sinusoids, with the same charac- 
teristics as the impulse response function of the original system. Thus, cross-correlation 
functions can be used as impulse response functions in time-domain modal parameter esti- 
mation schemes. This is more clearly seen after further simplification using a new constant 
multiplier (Gf), as follows: 


r = l m r i o d 


VK T 


to r d T + Qr 


James et al. [1995] present more definition of the intermediate steps in this derivation and the 
results of analytical verification checks of the derivation. 
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Comparison of NExT Test Data from a Rotating VA WT 
to Simulated Data 

An interesting comparison was performed between NExT data processing results and 
simulated data and using the known modal parameters. For this activity, a simulation code 
(VAWT-SDS), was used to compute the time domain response of the DOE/Sandia 34-m 
VAWT test bed (Figure 2-2) during rotation in turbulent wind [Dohrmann and Veers, 1 989]. 
The structural model (including rotational effects) used in VAWT-SDS was available to cal- 
culate the analytical modal frequencies and analytical modal damping. VAWT-SDS was 
used to generate analytical data, which were then input to NExT. The results were compared 
to the known frequencies and damping information to test the capabilities of NExT. 

Simulated data were generated for the 34-m testbed using a 30 rpm rotation rate and 20 
mph (8.9 m/s) turbulent winds with 15% turbulence intensity. Stiffness proportional damp- 
ing, sufficient to produce a damping ratio of 0.2% at 1.4 Hz was added to the model. Time 
histories of 2,048 points for eight strain gauge outputs were generated using a step size of 
0.04 s. Sensor noise was simulated by adding a white-noise signal to each simulated time 
history. The standard deviation of this additive signal was 2% of the standard deviation of 
each time history. The analytical modal frequencies and damping ratios were also calculated 
by extracting the complex eigenvalues from the structural matrices used in the VAWT-SDS 
code. VAWT-SDS used the Newmark-Beta numerical integration scheme, and the approxi- 
mations inherent in this procedure produced period elongations. The frequency shifts created 
by numerical integration were calculated and a correction was added to the analytical values 
[James et al., 1996]. 

NExT was used to estimate modal frequencies and damping ratios from the simulated 
data so that the NExT results and the analytical values could be compared. Table 11-2 shows 
the results of the comparison. With NExT, it was possible to correctly extract the frequen- 
cies and damping from the data with the exception of the first two modes at 1 .27 Hz and 1 .35 
Hz. These two modes are very closely spaced, making it difficult to obtain accurate frequen- 


Table 11-2. 

Comparison of NExT With Simulated Results [James et al., 1993] 


Mode 

Frequency (Hz) 

Damping (%) 

Simulated 

NExT 

Simulated 

NExT 

1st Flatwise Antisymmetric 

1.27 

1.31 

0.2 

0.4 

1st Flatwise Symmetric 

1.35 

1.32 

0.2 

0.3 

1st Blade Edgewise 

1.59 

1.59 

0.3 

0.3 

1st Tower In-Plane 

2.02 

2.01 

0.3 

0.4 

2nd Flatwise Symmetric 

2.43 

2.44 

0.4 

0.5 

2nd Flatwise Antisymmetric 

2.50 

2.50 

0.4 

0.4 

1st Tower Out-of-Plane 

2.80 

2.80 

0.3 

0.5 

2nd Rotor Twist 

3.39 

3.39 

0.5 

0.6 

2nd Tower In-Plane 

3.46 

3.45 

0.5 

0.4 

3rd Flatwise Antisymmetric 

3.65 

3.63 

0.5 

0.4 

3rd Flatwise Symmetric 

3.73 

3.73 

0.6 

0.4 

2nd Blade Edgewise 

3.88 

3.87 

0.5 

0.3 
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cy or damping results from NExT. It should also be noted that the higher modes (3.65 Hz, 
3.73 Hz, and 3.88 Hz) have NExT-estimated damping ratios that are lower than the specified 
damping ratios. The amplitudes of these modes are low compared to the noise level, which 
affected the estimates. This finding suggested the need for longer time histories to improve 
the accuracy. Also, this comparison shows the validity of natural wind excitation as a tool for 
wind turbine testing, as well as pointing out the issues (such as period elongations) associated 
with numerical integration schemes used in analysis codes. 

Comparison of NExT to Non- Rotating Step Relaxation Excitation 

Another insightful comparison is between NExT processing using wind excitation and 
traditional analysis using step-relaxation excitation. For this comparison, a FloWind Cor- 
poration 19-m VAWT in Altamont Pass, CA, was tested using step-relaxation modal testing 
techniques during quiescent daytime winds. NExT was then used during periods of more 
substantial nighttime winds (above 7 m/s or 16 mph). The turbine was parked (nonrotat- 
ing) during all testing. Accelerometers were used to measure the response at predetermined 
locations on the turbine. This allowed a comparison between modal parameters estimated 
by NExT and modal parameters estimated using conventional testing techniques (step- 
relaxation testing). 

Table 11-3 compares the modal frequencies and modal damping ratios of the 19-m 
VAWT as determined from conventional testing and from NExT. The two methods pro- 
duced estimates of the modal frequencies that are in good agreement, particularly in view of 
the temperature difference between day and night. The average difference for the ten modes 
is only 0.5%. Also, the modal damping ratios of all six of the tower modes (rotor twist, tower 
in-plane, tower out-of-plane) are very similar. However, the modal damping ratios of all 
four blade flatwise modes (flatwise symmetric and flatwise antisymmetric) are substantially 
higher from NExT estimates. 


Table 11-3. 


Comparison of NExT With Step-Relaxation Excitation [James etal., 1993] 


Mode 

Frequency (Hz) 

Damping (%) 

Step Relax 

NExT 

Step Relax 

NExT 

1st Rotor Twist 

2.37 

2.38 

0.2 

0.1 

1st Flatwise Antisymmetric 

2.48 

2.49 

0.2 

1.3 

1st Flatwise Symmetric 

2.51 

2.51 

0.1 

1.4 

1st Tower Out-of-Plane 

2.72 

2.76 

0.4 

0.4 

1st Tower In-Plane 

3.11 

3.15 

0.4 

0.4 

2nd Tower Out-of- Plane 

4.53 

4.53 

0.1 

0.1 

2nd Flatwise Antisymmetric 

5.30 

5.31 

0.3 

0.8 

2nd Flatwise Symmetric 

5.64 

5.65 

0.1 

0.6 

2nd Rotor Twist 

6.59 

6.62 

0.1 

0.1 

2nd Tower In-Plane 

6.64 

6.71 

0.3 

0.6 
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Subsequent analysis suggested that these differences were due to a drag phenomenon 
experienced as the largest cross-section of the blade oscillated normal to the stronger night- 
time winds [James et al., 1995]. This exercise showed the utility of performing modal test- 
ing in an operating environment to allow the extraction of the total damping (structural and 
aeroelastic). This test series also pointed out the reduced infrastructure associated with the 
natural excitation test while clarifying the utility of traditional testing to separate out the dif- 
ferent physical phenomena acting on wind turbine structures. 

Application of NExT to a Rotating VAWT 

The significant utility of wind excitation is the ability to extract structural dynamics in- 
formation from a rotating system. NExT was used to extract modal frequency and damping 
using data from the DOE/Sandia 34-m testbed during rotation at 0, 10, 15, 20, 28, 34, and 38 
rpm. The wind speed during these tests was also approximately 10 m/s (22 mph). Twelve 
strain gauges (accessed through slip rings) were used as the sensors in all of these analyses. 
Time history length, sample rate, correlation block size, and the time domain analysis pa- 
rameters were all the same as provided in the previous section. Figure 1 1-10 is a plot of the 
rotation-rate dependent modal frequencies of this turbine from analysis and from experiment 
using NExT. 


34 METER TEST BED FREQUENCIES 



Figure 11-10. Modal frequencies versus rotor rotation rate for the Sandia/DOE 34-m 
VAWT. [James etal., 1993] 
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Figure 11-11 shows a plot of modal damping ratios, as calculated with NExT, versus 
turbine rotation rate for the blade flatwise modes. The damping ratios of these modes gener- 
ally increase significantly with turbine rotation rate. For example, the first flatwise damping 
ratio increases from 2% to 7%. The notable exception is the second blade flatwise mode 
that drops between 15 rpm and 28 rpm. Such a drop in damping ratio is likely due to modal 
coupling to a more lightly damped mode, since modal coupling varies with rotation speed. 
This data validates the need to understand the variation in parameters in operating conditions 
including rotation rate effects, especially in variable- speed machines or during start-up/ 
shut-down. 



Figure 11-11. Modal damping versus rotor rotation rate for the Sandia/DOE 34-m 
VAWT. [James etal., 1993] 

Application of NExT to HAWT Systems 

Table 11-4 shows the results of comparison between impact excitation and natural ex- 
citation on a parked HAWT. In this case the machine was an Atlantic Orient Corporation 
(AOC) 15/50, which is a three-bladed turbine with a 15-m rotor mounted downwind of the 
tower. The results show a very close agreement between the measured frequencies and a 
typical agreement for the damping values [Rumsey et al., 1997]. 

A controlled yaw or upwind-rotor HAWT is the most direct application of NExT to a ro- 
tating turbine, as the kinematics of the system are primarily limited to rotation of the blades. 
This class of turbine shows the same dynamic characteristics as the rotating VAWT, including 
distinct natural frequencies, environment- dependent damping, and rotation-dependent forced 
harmonics. NExT has been applied to field results from a Northern Power Systems 100-kW 
HAWT, which has a two-bladed teetering rotor, 17.8 m in diameter [James, 1994], Figure 11- 
12 displays some representative data from this analysis. In this case the damping in the 5.34 Hz 
edgewise mode is plotted versus wind speed. The damping is seen to have a clear and increasing 
trend as the wind speed increases. 
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Table 11-4. 

Comparison of Two Modes from a Parked AOC 15/50 HAWT 
Using NExT and Impact Excitation [Rumsey et al., 1997] 


Mode 

Frequency (Hz) 

Damping (%) 


Impact 

NExT 

Impact 

NExT 

Teeter 

3.19 

3.18 

1.18 

1.58 

Umbrella 

3.72 

3.73 

1.11 

1.27 
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DAMPING VS. WIND SPEED FOR 5.34 Hz MODE 



WIND SPEED (MPH) 


Figure 11-12. Damping versus wind speed for the 5.34 Hz edgewise mode of the North- 
ern Power Systems 100-kW HAWT. [James, 1994] 


A more challenging dynamical system is a free-yaw HAWT with its rotor mounted down- 
wind of the tower. This turbine configuration can be simpler to build, but the uncontrolled 
yaw degree-of-freedom couples to the other in-plane dynamic modes and creates entire fami- 
lies of frequency-dependent phenomena. Figure 11-13 shows an example from the free-yaw 
26-m diameter Advanced Wind Turbines AWT-26, which is a downwind teetering-rotor ma- 
chine [Malcolm and James 1995]. The expected per-rev harmonic forcing function frequen- 
cies are plotted as the dotted lines. The solid lines denote the analytically-predicted natural 
frequency families. These frequencies are only distinct in the non-rotating condition. 

As the machine rotates, natural frequencies split into multiple rotation-rate dependent 
frequencies. NExT was used to process field data from this HAWT as it was rotating at 58 
rpm. There is a slight mismatch between the analysis and the field test data. This mismatch 
illustrates the utility of extracting structural dynamics data from an operating machine for 
comparisons to analytical data. 
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Figure 11-13. Periodic natural frequency families including NExT-derived field data 
vs. rotor speed for the AWT-26 HAWT. [Malcolm and James, 1995] 
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Classical Aeroelastic Flutter Analysis for HA WTs 


Background 

Flutter is a self-starting and potentially destructive vibration where aerodynamic forces 
on an object couple with a structure’s natural mode of vibration to produce rapid periodic 
motion . Flutter can occur in any object within a strong fluid flow, under the conditions that 
a positive feedback occurs between the structure’s natural vibration and the aerodynamic 
forces. That is, when the vibrational movement of an object like a flexible airfoil increases an 
aerodynamic load which in turn drives the object to move further. If the energy during the 
period of aerodynamic excitation is larger than the natural damping of the system, the level of 
vibration will increase. The vibration levels can thus build up and are only limited when the 
aerodynamic or mechanical damping of the object match the energy input. This often results 
in large amplitudes and can lead to rapid failure. 

Although classical aeroelastic flutter has generally not been a driving issue in utility- 
scale wind turbine design, one case in which classical flutter was observed involved a vertical- 
axis wind turbine (VAWT) turning in still air [Popelka 1982], The rotor was purposefully 
motored at ever-increasing speeds until the flutter boundary was breached and dramatic clas- 
sical flutter oscillations were observed. Flutter occurred at approximately twice the design 
operating speed of the rotor. 

For very large horizontal-axis wind turbines (HA WTs) fitted with relatively soft (flex- 
ible) blades, classical flutter becomes a more important design consideration. Innovative 
blade designs involving the use of aeroelastic tailoring, wherein the blade twists as it bends 
under the action of aerodynamic loads to shed loads resulting from wind turbulence, increase 
the blades proclivity for flutter. 

Analytical Model of Flutter 

The analysis of classical flutter in wind turbines necessitates the use of unsteady aerody- 
namics. As pointed out by Leishman [2002], for horizontal axis wind turbines there are two 
interconnected sources of unsteady aerodynamics. The first is a result of the trailing wake of 
the rotor and is addressed by investigating the interactions between the rotor motion and the 
inflow. The second, which will be the focus here, is due to the shed wake of the individual 
blades and can be addressed using techniques developed for analyzing flutter in fixed-wing 
aircraft. 

To simplify the analysis, the rotor is assumed to be turning in still air. Because there is 
no wind inflow, unsteady aerodynamics caused by the trailing wake can be neglected. Con- 
sequently, the aerodynamic behavior of a single blade is similar to that of a fixed wing with a 
free stream velocity that varies linearly from the root to the tip, assuming that the shed wake 
of the preceding blade dies out sufficiently fast so that the oncoming blade will encounter 
essentially still air on a hub fixed in space. Focusing on aeroelastic stability associated with 
the shed wake from an individual HAWT blade, the technique developed by Theodorsen 
[Theodorsen 1935, Fung 1969, Dowell 1995] for fixed wing aircraft has been adapted for use 
with HAWT blade flutter [Lobitz 2004], 

The Theodorsen technique specifically addresses classical flutter in an infinite (i.e. two 
dimensional with no end effects) airfoil undergoing oscillatory pitching and plunging motion 
in an incompressible flow, as illustrated in Figure 1 1-14. The airfoil’s pitching motion is rep- 
resented by the angle a and its plunging motion by the vertical translation h. L represents the 
lift force vector positioned at the quarter-chord, M is the pitching moment about the elastic 
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axis and U is the free-stream wind velocity. The origin of the coordinate system is positioned 
along the chordline at the section elastic axis. Note that a is defined as the fraction of b (the 
half-chord) that the elastic axis is aft of the mid-chord. Thus in Figure 1 1-14, a is negative 
since the elastic axis is ahead of the mid-chord. 



Figure 11-14. Schematic diagram of motions of a blade cross section assumed in the 
application of the Theodorsen technique for analyzing the flutter of a fixed wing and 
adapted for wind turbine blades. 


The blade is assumed to be simultaneously pitching and plunging in an oscillatory fash- 
ion at a modal natural frequency to, as follows: 

h= hoe l<ot and a = aoe‘ wt (11-55) Q5 


where 

h = vertical displacement (m) 
h Q = complex constant (m) 
to = modal natural frequency (rad/s) 
t = elapsed time (s) 
a = pitch angular displacement (rad) 
a 0 = complex constant (rad) 

According to the Theodorsen model, the lift force, L, and the pitching moment, M, are given 
by Equations (11-56), where a, b, d r and d 2 are as shown in Figure 11-14. 


( 


L= 27tp U 2 b\ 




io)C(k) , „ , 

— — — h 0 + C(k) (Xo + [1 + C(&)(1 - 2a)] — (Xo 

to 2 b to 2 b 2 a 

2U 2 h ° + 2 U 2 a ° 




J 


(11 -56a) 
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M = 2npU 2 b - 
where 

L = 
M = 
U = 
C(k) = 
k = 

The Theodorsen 

where H denote 
graphically in Fi; 

Chapter 11 

r \ 

, i(oC(k ) , n itob , io)b 

d\ jj ho + C(k) ao + [1 + C(k)( 1 Oq + d 2 a o 

00 2 ab 2 (\ 2 \ ® 2 ^ 3 

2Ui h0 + \S +a ) 2f/2 a ° J 

(ll-56b) 

lift force intensity; force per unit length of span (N/m) 

pitch moment intensity, moment per unit length of span (N-m/m) 

resultant airspeed (m/s) 

complex-valued Theodorsen function 

reduced frequency (rad) 

function, C(k) and the reduced frequency, k, are defined as follows: 

k=<nb/U (ll-57a) 

H^(k) 

C(k) = (2) 1 ( \ 2) (11 -57b) 

H\\k) + iH (2 \k) 

s the Hankel function. The real and imaginary parts of C are displayed 
gure 11-15. 
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Figure 11-15. Real and imaginary parts of the Theodorsen function C(k) where k is the 
reduced frequency. 
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In order to incorporate Equations (11-56) in a finite element procedure for subsequent 
complex eigenvalue analysis, they must be recast in a pseudo time domain form for develop- 
ing contributions to the finite element mass, stiffness and damping matrices. This can be 
accomplished by leaving the Theodorsen function as is, but using the explicit to’ s in the equa- 
tions to construct time derivatives, producing the following equations: 


L = 2npU 2 b ( Ci ^-h 


h + C(k ) a + [ 1 + C(k)( 1 - 2a)] j (H-58a) 


M = 2npU 2 b 


d x ^-h + C(k)a + [1 + C(k)(l -2a)]^<x + d 2 ~a 


ab 2 - 
+ 2 U^ h 


9 .. 

2U 2 a 


(11 -58b) 


Now, using Equations (11-58) with the principle of virtual work, contributions to the finite 
element stiffness, mass and damping matrices can be developed that include the complex- 
valued Theodorsen function. 

However, before that is done, it is noted that the HAWT blade does not conform to the 
configuration of an infinitely long uniform wing to which the above equations apply. Rather, 
such quantities as the semi-chord, b, and the resultant velocity, U, vary significantly along the 
span of the blade. Moreover, the lift curve slope, which in the above equations is assigned the 
theoretical value of 2 tc corresponding to a flat plate, varies with blade span. These variations 
can be approximated by assembling a conglomerate of uniform blades, wherein the above 
equations are assumed to be applicable incrementally. Specifically, the quantities mentioned 
above are represented by a linear variation over the length of the element and included within 
the integral over the element length associated with the principle of virtual work. 

Significant detail can be incoiporated in this analysis by refining some of the variable 
quantities mentioned above. For example, inboard stall for a twisted blade turning in still 
air might be accommodated by reducing the lift curve slope in that region of the blade. 
More complicated inflow variations can also be accommodated if they are known a priori. 
However, in comparison with other approximations that have been made in this analysis, 
these types of refinements are deemed excessive. The simple aeroelastic stability analysis 
presented here is meant to serve primarily as a common sense check for use during the blade 
design process. 

In matrix notation, the finite element equation to be used for investigating the aeroelastic 
stability of the wind turbine blade is as follows: 


[M + M a (Q)] {«} + [C c (Q) + C a (( o, Q)] {ii} 

+ [K(u 0 , Q) + K tc + K CS (U) + K a (to, Q)] {«} = 0 


(11-59) 


where M is the conventional mass matrix and A"(u 0 , Q) is the stiffness matrix, centrifugally 
stiffened commensurate with displacements, u g , resulting from centrifugal loads correspond- 
ing to the rotor rotational speed, Q. The displacement, u, and its time derivatives represent 
motion about this centrifugally loaded state. The a and h degrees of freedom of Equations 
(11-58) are included in u. 
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The blade is modeled with NASTRAN tapered beam (CBEAM) elements (developed 
by the MSC Software Corp.) that have six degrees of freedom at either end. These elements 
also have provisions for shear deformation and warping of the cross section, although these 
features were not used in this analysis. Additional terms are added to the various matrices 
to provide for rotating coordinate system effects. These are the Coriolis and centrifugal 
softening terms, C c (Q) and ^(Q), respectively, that will be discussed in detail later. The 
aeroelastic matrices, M a (Q), C a (co,Q) and ^(<n,Q), all depend on Q, since the free-stream 
velocity at any blade radius is governed by it (it is assumed that the rotor is turning in still air). 
Linear shape functions are used in the development of these aeroelastic matrices. 

As a demonstration of the method used to generate the aeroelastic matrices, a single 
beam element’s contribution to K a is presented in Equation (11-60). The element has a 
length, Z, and is bounded by nodes 1 and 2. The stiffness contribution is due to the a terms 
of Equations (11-58) and can be developed with the principle of virtual work using only a 
virtual displacements. Thus, for a varying linearly along the element (i.e. a(x) = a/l-xA) + 
a^/l ) the contribution of this single element is as follows: 




(1-f) 2 f(l-f) 

Li 

Mi 

i 

■ = p J ao(x)U 1 (x)b(x)C(k(x)) 
0 

f(l-f) (f) 2 

*(*) (1 - f) 2 (1 " 7 ) 

m 2j 


x*)f( 1-f) <fi(*)(f) 2 


OC] 

a.2 


(11-60) 


Here the quantity 2n of Equations (1 1-58), which represents the lift curve slope for a flat 
plate, is replaced by a Q (x) representing the lift curve slope of an airfoil. This slope varies 
along the length of the element as does the free stream velocity, U, the semichord, h, the 
Theodorsen function C, and the distance from the elastic axis to the quarter-chord, d y The 
integral is evaluated numerically. Elemental contributions to M a and C a are developed in an 
entirely similar manner. 

Commensurate with the use of the NASTRAN CBEAM element above, the NASTRAN 
commercial finite element software is used for this aeroelastic stability investigation. The 
contributions to the stiffness, mass, and damping matrices discussed above are supplied by 
means of a NASTRAN input option. NASTRAN can accommodate non-symmetric, com- 
plex-valued matrices as are required in this effort, and it provides a number of complex 
eigenvalue solvers for the stability analysis. 

In addition to £2, the aeroelastic matrices, C o (go,Q) and K a (oa,Q), also depend on oo, the 
natural frequency of the mode shape of interest, which occurs in the argument of the The- 
odorsen function. Since this frequency is unknown at the onset of the computations an itera- 
tive process is required for obtaining accurate results. The iterative procedure developed for 
the aeroelastic stability analysis of the rotor blade is composed of the following steps: 

1. Select a value for Q . 

2. In a quasi-static NASTRAN run, create K(u 0 ,Q) for subsequent eigenvalue 
analysis. 
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3. Provide an initial guess for co or update it from the prior calculation. 

4. Using a NASTRAN complex eigenvalue solution procedure, compute modes, fre- 
quencies and damping coefficients. 

5. Select a mode with a small or negative damping coefficient and return to step 3 with 
corresponding frequency update. 

6. When the prior updated frequency is sufficiently close to the subsequently com- 
puted one, either suspend computations or modify Q and return to step 1. 

The end goal in this process is to identify the eigenmode that exhibits a negative damp- 
ing coefficient for the lowest rotor rotational speed. This speed is designated as the classical 
flutter speed for the blade. A comprehensive approach would be to sequentially investigate 
the stability of each mode by choosing its frequency as the initial guess for to, and tracking 
that frequency over a range of increasing Q. The flutter speed would then correspond to the 
lowest value of Q where the damping becomes negative over the set of modes investigated. 

In practice, however, the first torsional blade mode produces the lowest values of Q, coupling 
with other blade modes (primarily flapwise modes) as the rotational speed increases. Of 
course it is always prudent to check other modes that appear suspicious during the course of 
the analysis. 

Sample Flutter Calculation 

After developing some confidence in the above computational procedure by comparing 
its results to those of Theodorsen, it was applied to the large relatively soft wind turbine blade 
shown in the wire-frame illustration of Figure 1 1-16. This blade, termed the baseline blade, 
was designed as part of the WindPACT Rotor Design Study [Malcolm and Hansen 2002]. 

Some characteristics of this blade and the associated rotor are listed in Table 1 1 -5. Q6 



Figure 11-16. Wire-frame illustration of the 1.5 Mw baseline WindPACT blade. 

[Malcom and Hansen 2002] 
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Table 11-5. 


Characteristics of the Baseline 1.5 MW WindPACT Blade and Rotor 


Characteristic 

Units 

Value 

Rated power 

MW 

1.5 

Rotor diameter 

m 

70 

Max rotor speed 

rpm 

20.5 (0.342 Hz) 

Max tip speed 

m/s 

75 

Blade coning 

deg 

0 

Max blade chord 

m 

8% of radius 

Radius to blade root 

m 

5% of radius 

Rotor solidity 


0.05 

Blade mass 

kg 

4,230 

Hub mass 

kg 

15,104 

Total rotor mass 

kg 

32,016 

1 st flap wise frequency 

Hz 

1.233 (1.199) 

1 st edgewise frequency 

Hz 

1.861 (1.714) 

2 nd flapwise frequency 

Hz 

3.650 (3.596) 

1 st torsional frequency 

Hz 

9.289 (9.846) 


Although not readily apparent from Figure 11-16, the blade has a modest amount of 
twist, from 11.1 deg. at the root to 0.0 deg. at the tip. The modal frequencies in parentheses 
in this table were computed with the software developed herein with Q set to zero and the 
Theodorsen function, C(k), set to 1 .0. They differ from the companion frequencies taken 
from the WindPACT study, which were computed using the ADAMS software (developed 
by the MSC Software Corp.). The two models are inherently different, yet the frequency dif- 
ferences are small. Only minimal effort was made to minimize these differences. 

In Figure 1 1-17 a planform view of the blade is presented showing the positions of the 
elastic axis and the axis representing the locus of the cross-sectional center of gravity. As 
the elastic axis is forward of the mid-chord, the quantity a in Equations (1 1-59) is negative, 
varying along the span from - 0.354 to - 0.284. Since in this implementation a constant value 
of a is required, a mid value of - 0.32 was selected. 


<D 


I 

£ 

u 



Distance from rotor center (m) 


Figure 11-17. Diagram of the planform of the baseline 1.5 MW WindPACT blade show- 
ing the locations of the elastic axis and the cross-sectional centers of gravity. 
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Flutter Analysis Results 

Aeroelastic stability predictions with the rotor turning in still air were made for this base- 
line model. Results indicate that the onset of flutter occurs at a rotational speed of 43.4 rpm 
(0.723 Hz) which is 2.14 times the maximum design operating speed of the rotor. As shown 
in Figure 11-18, the mode shape associated with the onset of flutter contains significant 
amounts of edgewise (motion in the direction of the chord), flapping, and pitching motion. 
Its frequency of oscillation is 6.234 Hz, which is significantly lower than the 1 st torsional 
frequency of 9.289 Hz. 



Distance from rotor center (m) 
Amplitude 


Distance from rotor center (m) 
Phase (deg) 


Figure 11-18. Mode shape at the onset of flutter for the 1.5 MW baseline WindPACT 
uncoupled blade. 


In classical flutter, motion in the edgewise direction is generally minimal, and the large 
amount predicted here is thought to be due to the fact that the blade is twisted, and that the 
frequency of the second edgewise mode (6.153 Hz) is close to the flutter frequency of 6.234 
Hz. The amount of twist in the blade appears to have an important influence on the flutter 
mode shape. In fact, when flutter predictions are made for an entirely similar blade that is not 
twisted, edgewise motion in the associated flutter mode is minimal. 
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According to this analysis, the 2 nd flapwise mode combines with the 1 st torsional mode 
to form the flutter mode. This combination of modes was also observed by Hansen [2004] 
in predicting classical flutter in HAWTs with a software capability he has developed. In 
contrast, for fixed wing aircraft, it is the 1 st flapwise mode that combines with the 1 st torsional 
mode. 

Experimental Verification 

Unfortunately, no experimental data exist to corroborate these results. However, the 
flutter analysis technique described herein has been verified by comparison to the observed 
flutter in a VAWT [Lobitz and Ashwill 1985]. The experimental flutter frequency was mea- 
sured at 745 rpm [Popelka 1982], reasonably close to the predicted frequency of 680 rpm. 
The inclusion of a small amount of structural damping in the model greatly improved the 
agreement. The flutter mode shape was also well predicted using this technique. 
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